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Abstract 

Consider a two-dimensional continuous-time dynamical system, with an attracting fixed 
point S. If the deterministic dynamics are perturbed by white noise (random perturbations) of 
strength e, the system state will eventually leave the domain of attraction D.oi S. We analyse 
the case when, as e ^ 0, the exit location on the boundary dD. is increasingly concentrated near 
a saddle point H of the deterministic dynamics. We show using formal methods that the asymp- 
totic form of the exit location distribution on dD. is generically non-Gaussian and asymmetric, 
and classify the possible limiting distributions. A key role is played by a parameter //, equal to 
the ratio \\g(H)\/\^(H) of the stable and unstable eigenvalues of the linearized deterministic 
flow at if . If // < 1 then the exit location distribution is generically asymptotic as e ^ to a 
WeibuU distribution with shape parameter 2///, on the 0(e''''^) lengthscale near if. If// > litis 
generically asymptotic to a distribution on the 0(e^^^) lengthscale, whose moments we com- 
pute. The asymmetry of the asymptotic exit location distribution is attributable to the generic 
presence of a 'classically forbidden' region: a wedge-shaped subset of Q. with H as vertex, 
which is reached from 5, in the e ^ limit, only via 'bent' (non-smooth) fluctuational paths 
that first pass through the vicinity of H. We show that as a consequence the Wentzell-Freidlin 
quasipotential function W, which governs the frequency of fluctuations to the vicinity of any 
point a; in il and is the solution of a Hamilton-Jacobi equation, generically fails to be twice 
differentiable at a; = ii. This nondifferentiability implies that the classical Eyring formula 
for the small-e exponential asymptotics of the mean first exit time, which includes a prefactor 
involving the Hessian of W^ at a; = ii, is generically inapplicable. Our treatment employs 
both matched asymptotic expansions and probabilistic analysis. Besides relating our results 
to the work of others on the stochastic exit problem, we comment on their implication for the 
two-dimensional analogue of Ackerberg-O'Malley resonance. 



1 Introduction 

We consider the problem of noise- activated escape from a connected planar domain il C M^ with 
smooth boundary, in the limit of weak noise. If b = (b'), i = l,2isa smooth vector field on a 
neighborhood of the closure Q, we define the random process x^(t) = (xl(t)), i = 1, 2 by the Ito 
stochastic differential equation 

dxlit) = b\xM) + e'/' E <^'a(xM) dwjt) (1) 

a 

and an appropriate initial condition. Here Wa(t), a = 1,2, are independent Wiener processes and 
(^ = (cr'a) is a 2-by-2 noise matrix, like b a function of position x = (x'), i = 1,2. The associated 
diffusion tensor D = (D'^) is defined by 

D'^ =J2a\a^^, (2) 

a 

i.e., D = (Tcr^ . We assume strict ellipticity, i.e., that D is nonsingular on Q. and its boundary, and 
that its a; -dependence is smooth on a neighborhood ofCl. e > is a noise strength parameter. The 
subscripts on x\ and x,^ emphasize the e-dependence of the random process, which may be viewed 
as a dynamical system stochastically perturbed by noise. 

Of interest in applications is the case when Q. contains only a single stable fixed point S of 
the drift field b, and S serves as an attractor for the whole of Q.. If Q. is the entire domain of 
attraction of S, the boundary dQ. will not be attracted to S: it will be a separatrix between domains 
of attraction. This 'characteristic boundary case' is particularly important and difficult to study. 
We shall analyse this case, assuming that dO. is a smooth characteristic curve of b containing 
fixed points (alternating saddle points and unstable fixed points, as in Figure 1). We allow Q. to 
be unbounded. If the initial condition is aj^CO) = S and r^ is the first passage time from S to 
the boundary (i.e., the first exit time), we shall study the behavior of the exit location distribution 
Peix) dx = ? {x,XT() ^ X + dx] on dQ. in the e ^ limit, and the small-e asymptotics of the 
mean first passage time (MFPT) Er^. 

Exit problems of this sort, in more than one dimension, have a long history [65]. They arose 
originally in chemical physics [2, 4, 9, 29, 40], but occur in other fields of physics [64, 75] 
as well as in systems engineering [53, 76] and theoretical ecology [49, 57, 58]. In recent years 
two different approaches have been used: rigorous large deviations theory [14, 15, 26, 27] and 




Figure 1: The structure of the drift field b on Q., when Q. is bounded. The boundary dO., being 
a separatrix, is not attracted to the stable fixed point S, though all points in Q. are attracted to S. 
Saddle points (H) and unstable fixed points (U) alternate around the boundary. 



formal but systematic asymptotic expansions [49, 61, 65, 74]. The rigorous approach yields 
comparatively weak but still very useful results. In particular much light is thrown on exit 
problems by the Wentzell-Freidlin quasipotential [27], or classical action function, W : Cl ^ M"^. 
W(x) is best thought of as a measure of how difficult it is for the process x^(t) to reach the point x ; 
as e ^ the frequency of excursions to the vicinity of x will be suppressed exponentially, with rate 
constant W(x). As we review in Section 4, W(x) has an interpretation in terms of certain 'classical' 
trajectories extending from S* to a;, interpreted as the most probable (as e ^ 0) fluctuational paths 
from S* to a; [49,56]. 

Normally one expects that W will attain a minimum on dD. at one of the fixed points, in particular 
at a saddle point. This can be shown to imply that as e ^ 0, the exit location on dQ. converges 
in probability to the saddle point. Actually the behavior of W on dQ, and its consequences, are not 
yet fully understood. (For some partial results on the smoothness of W, see Day and Darden [19] 
and Day [18].) Recent treatments [16, 54] indicate that as e ^ it is possible for the exit location 
to converge to an unstable fixed point on dQ., due to local constancy of W on the boundary. In this 
paper we consider only models displaying the conventional behavior: as e ^ the exit location 
should converge to some distinguished saddle point H on dQ, due to W on dQ having a unique 
global minimum at H. 

The method of formal asymptotic expansions has provided evidence for an unusual phe- 
nomenon: skewing of the exit location distribution in the vicinity of the saddle point H. What is 



meant by skewing is that as e ^ 0, the exit location converges to H but its distribution p^ix) dx 
on d^ is not asymptotic to a Gaussian centered on H. (Skewing was discovered by Bobrovsky 
and Schuss [7]; for recent rigorous work, see Bobrovsky and Zeitouni [8] and Day [17].) This 
phenomenon has until now remained unclarified, though examples of skewed asymptotic exit lo- 
cation distributions have been given by Maier and Stein [55]. Using formal methods we show in 
this paper that skewing is a generic phenomenon. We also derive a general result, analogous to the 
central limit theorem, which characterizes skewing: A* e ^ 0, in any generic stochastic exit model 
(with characteristic boundary) of the above sort, the exit location distribution, on an appropriate 
e-dependent lengthscale near H, will be asymptotic to a non-Gaussian distribution that belongs 
to one of two well-defined classes. Which asymptotic distribution occurs is largely determined by 
the behavior of the model near H (i.e., the ratio // = \Xs(H)\/Xu(H) of the stable and unstable 
eigenvalues of the linearization of b at H, which we assume to be nonsingular, and to a lesser 
extent the value of D(H)). One of the two classes is the class of WeibuU distributions, which is 
familiar from statistics [3]. The asymptotic exit location distributions in the second class are more 
complicated, and we do not derive explicit expressions for them here. We do however provide 
an algorithm for computing their moments, in terms of the correlation functions of a conditioned 
three-dimensional Bessel process. 

Much previous work, in particular on physical applications [32], has dealt with a special case: 
when the mean drift 6' equals — D'^9j4>, for some smooth potential function $. (Here d^ = d/dx^; 
summation over repeated Roman indices is assumed henceforth.) In this case the point H will 
simply be a saddle point of $. This gradient (or 'conservative') case finds many applications in 
physics, but from a mathematical point of view is highly nongeneric. W can be solved for exactly 
(it equals 2$), and one can show formally that in the limit of weak noise the distribution of the 
exit location on dO. is asymptotic to a Gaussian centered on H, with 0(e^^^) standard deviation. 
It turns out that generic, nongradient drift fields, which find application in nonequilibrium statistical 
mechanics and stochastic modelling generally, are qualitatively different. 

In the light of the large deviations approach, the fact that generic drift fields give rise to 
non-Gaussian asymptotic exit location distributions is surprising. One might expect that on dO., 
W in general attains a quadratic minimum at H. If d^W/ds^(s = 0) equals cr"^ (s being the arc 
length along dQ., measured from H), the abovementioned exponential suppression as e ^ would 
presumably give rise to a factor exp(—s^/2(7^e) in the exit location density Pe(s). This is precisely 
a Gaussian centered on H, on the (9(e^/^) lengthscale. The generic case would therefore seem to 
resemble the gradient case. 



It is tempting to ascribe the discrepancy between this prediction and our results to what 
theoretical physicists would call a 'prefactor effect', i.e., the presence of subdominant (as e ^ 0) 
terms in the exit location density that are not included in the usual exponential factor arising from 
large deviations theory. Surprisingly, this is not the case. The prediction assumed that if W has a 
minimum at H, it has a quadratic minimum there. We shall see that generically, W is not even twice 
dijferentiable at H\ strikingly, d'^Wjds^ is generically discontinuous at s = 0. This discontinuity 
causes the exit location density Pe{s), even if (in the small- e limit) it is localized on the (9(e^/^) 
lengthscale near H, to have different Gaussian falloff rates as sje^l'^ -^ +00 and sje^l'^ -^ —00. 
This is the true cause of skewing. The abovementioned limiting WeibuU distributions, which are 
one-sided, arise when either d^W/ds^(0+) or d^W/ds^(0—) equals zero (a generic occurrence, 
when /i < 1); on account of the corresponding zero Gaussian falloff rate, they are localized on a 
larger lengthscale than (9(e^/^). The appropriate lengthscale turns out to be (9(e^/^). 

The cause of the generic discontinuity in d^W/ds^ at s = is explained in Section 5; it has 
an easy (indeed pictorial) interpretation. The reader may wish to glance ahead at Figure 3, which 
displays the typical behavior of the piecewise classical trajectories giving rise to W(x), for x in 
the vicinity of H. The difference between trajectories incident on boundary points with s > and 
those incident on boundary points with s < is apparent. There is a wedge-shaped 'classically 
forbidden' region emanating from H, which is reached only by piecewise classical (rather than 
classical) trajectories. This region includes all the boundary points to one side of H. 

In Sections 2, 3 and 4 we explain the aspects of our approach that make contact with previous 
treatments: the reduction of the exit location problem to the study of the principal eigenfunction of 
the forward Kolmogorov equation, and the approximation of this eigenfunction by a characteristic 
(ray) expansion [49]. The rays employed will be the classical trajectories of Wentzell-Freidlin 
theory, and the same action function W will appear. We obtain a system of ordinary differential 
equations which must be integrated along the rays [49, 55, 56] . Our equations extend those obtained 
by Matkowsky, Schuss and Tier [60, 61]; being ordinary rather than partial differential equations, 
they are more suited to numerical computation. A key role is played by a Riccati equation for the 
matrix of second derivatives Z = (didjW), i,j = 1,2. We comment on the fact that our equations 
are covariant: they transform systematically under changes of coordinates, and have a geometric 
(coordinate-free) interpretation. 

In Section 6 we consider the nongeneric case when (d^W/ds^)(0+) and (d^W/ds^)(0—) are 
equal, which gives rise to a Gaussian asymptotic exit location distribution. In Sections 7 and 8 we 
turn to the generic case, and analyse models with // < 1 and // > 1. We study the fine behavior 



of the process x^(t) and the function W in the vicinity of H, and by a combination of matched 
asymptotic expansions (involving the construction of a novel 'inner' approximation in a boundary 
layer near dO.) and probabilistic analysis, derive the possible asymptotic exit location distributions 
on the appropriate lengthscales. We also extend earlier work by showing how our system of 
coupled ordinary differential equations may be employed to compute the small- e asymptotics of 
the MFPT Er^. It is known that Er^ grows exponentially as e ^ 0, with rate constant W(H). 
But in the case of a characteristic boundary, the numerical computation of the pre-exponential 
factor in the MFPT asymptotics is more complicated than previous treatments [60, 61, 65, 69, 74] 
have suggested. Until now it has been generally believed that the pre-exponential factor may be 
computed from the Hessian matrix of W at the saddle point H. But generically this matrix of 
second derivatives does not exist! We show how this is due to the ray expansion becoming singular 
at H, and how previous treatments may be modified to take this phenomenon into account. 

In Section 9 we comment on the implications of our results for a related area: the solution 
of multidimensional singularly perturbed elliptic boundary value problems. In one-dimensional 
singularly perturbed boundary value problems with turning points, the phenomenon of Ackerberg- 
O'Malley resonance [1, 66] is well understood [42, 43], but the corresponding multidimensional 
phenomenon has been little explored. We are able to show that at least one case of two-dimensional 
Ackerberg-O'Malley resonance, which can be reduced to a stochastic exit problem on Q. (dQ. being 
characteristic, and containing fixed points) has unusual features. On account of skewing near the 
dominant fixed point and an anomalous exit location lengthscale, a singular perturbation expansion 
for the matching 'outer' solution in the body of Q. must in general employ irrational powers of the 
perturbation strength e. This is strikingly different from the case of a non-characteristic boundary, 
where an expansion in powers of e suffices [60, 61]. It also differs from the corresponding one- 
dimensional problem (on an interval whose endpoints are turning points), where subdominant terms 
in the outer expansion are exponentially suppressed. Our results indicate that multidimensional 
Ackerberg-O'Malley resonance is likely to be more complicated than one-dimensional resonance. 
Our overall conclusions appear in Section 10. 

2 Mathematical Preliminaries 

The random process x^(t), t > has generator 

£, = -(e/2)D'^d,d, - Vd, (3) 



with formal adjoint £* defined by 

Clp = -(e/2)d,d,[D'^ p] + ddb'pl (4) 

The density pix^t) of the probability distribution p{x,t)dx = P {x^(t) G x + dx} will satisfy 
p = —jC*p, the forward Kolmogorov (or Fokker-Planck) equation. So the spectral theory of 
the operators C^ and C* is of interest. If they are equipped with Dirichlet (absorbing) boundary 
conditions on d^ and Q. is bounded, it follows by standard methods [39] that they have pure point 
spectrum, with smooth eigenfunctions, and that the principal eigenvalues of C^ and C* (the ones 
with minimum real part) are real. Also, the corresponding principal eigenfunctions may be taken 
to be real and positive on Q.. If Q. is unbounded we assume these properties continue to hold; 
they will hold if the stochastic model is sufficiently stable at infinity. We write Xf\ A^'\ Xf\ . . . 
for the eigenvalues of C^, ordered by increasing real part, and uJi*, m^\ u^, . . . for the corresponding 
eigenfunctions. Since Xf^ is real, the eigenvalues of C* will be A^*'^ A^^^*, A^*, . . . ; we denote the 
corresponding eigenfunctions by u*', u^\ u^, . . . If (•, •) denotes the usual L^ inner product on Q., 
we choose the u" and u" to satisfy (u™, u") = 6.mn- 

It is noteworthy that if b' = —D'^dj^ for some function 4> on Q., then the operator C^ is related 
by a similarity transformation, known as the Liouville transformation, to a formally self-adjoint 
operator; equivalently, C^ has a domain of self-adjointness. Consequently in the gradient case, 
the eigenvalues A^"^ will be real for all e > 0. The Liouville transformation is widely used in 
the physics literature [10, 63]. But although the gradient case occurs frequently in applications 
(for example, in the modelling of overdamped [inertialess] motion in a conservative force field, 
by means of a Smoluchowski equation), it is not generic. 

As e ^ the operators C^ and C* degenerate to first-order operators Co and Cq, so analysis 
of the e ^ limit is a singular perturbation problem. Formally one expects that A^"^ -^ Aq"^ and 
u" -^ Uq, where Aq"^ and u^ satisfy the first-order differential equation Cqu^ = Aq^^Mq. It will 
not in general be possible for the function Uq to satisfy the Dirichlet boundary conditions on dO., 
so a boundary layer near dO. of thickness tending to zero as e ^ is expected to be present. The 
convergence u" -^ Uq will be uniform only on compact subsets of ^. 

The formal limiting eigenvalues Aq"^ are easily worked out. Let B = (B'j) be the linearization 
of the drift field b' at x = S, i.e., B'j = djb\S)\ consider the case when B is negative definite. 



Approximation of jC^Uq = Aq^^Mq near 5*, using b'(x) sa B'j(x^ — S^), gives 

- B',(x^ - S')d,u^, = X^^^ul (5) 

In this linear approximation the uq will be homogeneous polynomials in the normalized coordinates 
x' — S'; Uq will be constant and the u^, n > 1, will be of higher degree. For each choice of 
positive degree the system (5) may be solved for the polynomial approximate eigenfunctions and 
their eigenvalues. In the gradient case, for example, the eigenvalues turn out to be of the form 
ni |6*^'^| + n2\b^^^\, rii G N, for b'^^\ b^^^ the two (negative) eigenvalues of B. 

The fact that Uq is always a constant function obviously holds beyond the linear approximation; 
the corresponding limiting eigenvalue \q^ equals zero. But the principal eigenfunction u*' of jC* is 
less well-behaved than u*^' as e ^ 0. An analysis similar to the above indicates that the u", for all n, 
are localized on the (9(e^/^) lengthscale near S. Their limits as e ^ are accordingly singular. 
The asymptotics of the principal eigenfunction v^ will be studied in detail in Sections 3, 6, and 7. 

It is easier to justify rigorously this formal analysis, and control the e ^ limit of Xf^ and u", 
when the boundary of Q. is non-characteristic: the integral curves of b cross dO. at all points, so that 
do. as well as Q. is attracted to S. In this case Friedman [28] proved that Xf^ -^ exponentially, 
and bounded it above and below; also Devinatz and Friedman [21] proved that u^ converges to a 
constant uniformly on compact subsets of Q.. In the one-dimensional case asymptotics of higher 
order can be obtained; see de Groen [20]. The behavior of the higher (n > 1) eigenmodes seems 
not to have been rigorously investigated, but exponential convergence of A^"^ to Aq"^ is expected 
as e ^ 0. 

In the case of a characteristic boundary there are some rigorous results of Eizenberg and 
Kifer [26] on the e ^ asymptotics of Xf\ but they do not apply when the boundary contains 
fixed points. Most rigorous work on the characteristic boundary case is probabilistic rather than 
analytic, so we remind the reader of the connections between the operators Ce,C* and the stochastic 
exit problem. Since the forward Kolmogorov equation p = —C*p may be written as a continuity 
equation p + diJ' = 0, with the probability current 

J' = -(e/2)d,[D'^p] + b'p, (6) 

the influx of probability into dO. (on which p = on account of absorption) has density 
—(€/2)dj[D^^ p](x) at point x G dO.. Accordingly the exit location density p^, as a function 



of a; G do., will satisfy 

p,(x) = j^ J\x,t)n,(x) dt = -{ell)n,{x)d, \^D'^ j^ p{-,t) dt^ (x) (7) 

where Uiix) is the outward normal to d^ dXx. In the special case when x,XO) = y this yields 

P.,y(^) = -(e/2)n,(x)d, \^D'^ f^^[cxp(-C:t)6y] dt^ (x) 

= -(e/2)n,(x)d,{D'^[(C:r'Sy]}(x) (8) 

the second subscript on p^ now denoting the initial condition. Here 6y is a Dirac delta function. 
But we have 

oo , 



/ J \ e 
n=0 



since the right-hand side is the integral kernel of (£*) ^ Substitution into (8) yields 

oo . 

V.^,y{x) = -(e/2)n,(x) V (A^"^)*) " <(y)a,[i)'-'<](a;). (10) 



n=Q 



Equation (10) is the point of contact between the analytic and probabilistic approaches to the 
stochastic exit problem. In the analytic approach asymptotic expansions, either formal or rigorously 
justified, are constructed for the eigenvalues and eigenfunctions. In the probabilistic approach the 
e ^ asymptotics of the exit location measure Pe,y{x) dx and the corresponding MFPT Ej^r^ are 
studied directly. 

The rigorous probabilistic approach was pioneered by Wentzell and Freidlin; in the case of non- 
characteristic boundary their results are as follows [27]. If the classical action function W attains a 
global minimum on d^ at some point x* , for any y ^ Q. the first hitting location on dQ. is as e ^ 
increasingly concentrated near x* , i.e., p^^y -^ 6x* as e ^ 0. Moreover EyT^ ~ exp(W(x*)/e), 
e ^ 0; the action W(x*) is interpreted as the asymptotic exponential growth rate of the MFPT. 
These results have been extended to the case of a characteristic boundary d^ by Day [15]. 

The lack of ^/-dependence in the e ^ limit has a probabilistic interpretation: irrespective of 
the choice of starting point y ^ Q., for sufficiently small e a typical sample path of the process 
x^t) tends to flow toward the attracting point S, and remain near S for a substantial (exponentially 
long) time, before experiencing a fluctuation large enough to drive it to dQ.. (With high probability 
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this will be along a Wentzell-Freidlin classical trajectory.) The lack of ^/-dependence also has an 
interpretation in terms of the series expansion (10). Since A^"^ decays exponentially as e ^ and 
the other A^"^ are expected to converge to nonzero values, up to exponentially small relative errors 
it should suffice to keep only the n = term in (10). That is, 

Pe,y(^) ~ -(e/2) (Af )"' u%)n,(x)d,Wv^](x), e ^ 0. (11) 



The asymptotic constancy of u'^ on Q., from this point of view, is the cause of the asymptotic 
^/-independence of the measure Pe,yix) dx on dQ.. Alternatively the boundary layer displayed by 
M° near d^ may be viewed as a consequence of the fact that when the starting point y is sufficiently 
near dQ., exit from Q. is likely to occur immediately rather than after a long sojourn in the vicinity 
ofS. 

The asymptotic validity of the reduction from (10) to (1 1) has not been rigorously proved, but 
is assumed (often implicitly) in most formal treatments of the exit problem. A related assumption, 
which we shall also make, is that up to exponentially small relative errors 

E.T, ~(Af)"\ e^Q (12) 



for any y ^ Q, the subscript y on the expectation denoting the initial condition for the process. 
This is similar to assuming that as e ^ the integral kernel of the evolution operator exp(— £*t), 
which equals the infinite sum 



^exp(-A('^)*t) u:(y)v:(x), (13) 

n=0 

may be approximated to high accuracy at large t by the n = term. If so the distribution of r^ for 
small e will be essentially exponential, with expectation ~ ( Af M as in (12). The assumption that 



only the n = term in ( 1 3) is significant at large t is justified by the fact that if the higher eigenvalues 
X^^\ n > I tend to nonzero limits as e ^ 0, the n > 1 terms, considered individually, should unlike 
the n = term decay on an 0(1) timescale rather than an exponentially long timescale. As e ^ 
there occurs a separation of timescales . 

The assumption (12) ties asymptotic expansions and Wentzell-Freidlin theory firmly together: 
it permits the identification of W(x*) with the asymptotic exponential decay rate of A^^ as e ^ 0. 
Note that whether or not this assumption is made, for any e > the eigenvalue A^"^ always has a 
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rigorous interpretation as the exponential decay rate of P {r^ > t} ast ^ oo [39] . Similarly for any 
e > 0, the eigenfunction u*' of C* corresponding to A^"^ has an interpretation as a quasistationary 
density. As t ^ oo, the n = term in (13) will dominate, and irrespective of the choice of initial 
condition /o(-,0), 

/>(-,t) = exp(-£:tV(-,0) (14) 

will be asymptotically proportional to uj^. 

At this point we introduce the singularly perturbed boundary value problem mentioned in 
Section 1 . Consider the problem 

Cu, = in Q., u, = f on dQ., (15) 

in which / : d^ ^ IR is a specified smooth function. It is easily checked that by Green's Theorem, 
for any e > the problem has solution 

u,{y) = -(e/2) / f{x)d,{W''{C:)-'8y]}{x)n,{x)dx. (16) 

That is, by (8), 

Ue(y)= f(x)p,^y(x)dx = Eyf(x,(T,)). (17) 

Jxedo. 

In other words u^(y) equals the expected value of / at the first hitting point x of the random process 
x^(t) on do., provided that the starting point ^^(O) equals y. This is an analytic rederivation of a 
well-known fact from the theory of random processes. In particular if p^y -^ 8x*, as will be the 
case if W attains a minimum on dQ. at x* , u, on Q. will tend as e ^ to a constant, namely fix*). 
And by (1 1) and (17), it will be possible to work out the small-e asymptotics of u, on Q. from the 
small-e asymptotics of u*' near dQ.. 

There are a handful of rigorous results on the boundary value problem which suggest that 
the asymptotic validity of the reduction from (10) to (11) is correct. Kamin [38] showed in the 
case of non-characteristic boundary that even if W does not attain a unique minimum on dQ., the 
solution u, 'levels' as e ^ 0, i.e., tends to a constant. Eizenberg [25] has shown that this levelling 
occurs exponentially rapidly, with rate constant that may be taken arbitrarily close to the infimum 
of W on dQ. This is consistent with the n > 1 terms in (10) being smaller than the n = 
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term by an exponentially small O(Xf^) factor. It is also known that u^ will have a boundary layer 
near dO., of width 0(e) if the boundary is non-characteristic and of width (9(e^/^) if the boundary is 
characteristic [33-35, 37]. This boundary layer may be viewed as arising from the boundary layer 
of the dominant eigenfunction u^. The boundary value problem will be studied further in Section 9. 

3 Matched Asymptotic Expansions 

It is clear from the last section that the weak-noise (e -^ 0) asymptotics of the stochastic exit 
problem are determined by the quasistationary density, i.e., the principal eigenfunction u*' of the 
forward Kolmogorov operator C*. The principal eigenfunction uj? of the backward operator C^ is 
expected to level exponentially rapidly as e ^ 0, so for all starting points y ^ O. the expression (11) 
for the e ^ asymptotics of the density p,^y of the exit location measure on 9Q may be written 

p.^yix) oc n,ix)d,Wv%x), e -^ 0, (18) 

where ni(x) is the outward normal to dQ. at x. The constant of proportionality here, which is 
independent of the starting point y, is fixed by the normalization condition 

/ p,Jx)dx = l. (19) 

Jxeo. 

We shall take y = S henceforth, though any other y ^ Q. could be chosen, and shall often drop the 
subscript. The interpretation of the asymptotic j/ -independence is as discussed in the last section. 
Pe^s is by (18) asymptotic to the absorption location density of the eigenmode v'^ on d^; also, 
the corresponding eigenvalue Af ^ of C* may be viewed as the absorption rate of the mode. If Jq is 
the probability current arising from Vq, i.e., 

Jl, = -(e/2)d,Wv'^] + kv'^ (20) 

then 

Af = / J^(x)n,(x)dx (21) 

Jxedo. 

provided that u*' is normalized to total unit mass. If this normalization condition does not hold. 
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(21) must be replaced by 



1(0) _ JxedQ. 'Jo\X)'>T-i\X) dx 
Ixen ^%x) dx 



-1 



Since we are assuming that (up to exponentially small relative errors) Est^ ~ ( A^"M , the MFPT 
asymptotics are determined by the asymptotics of v^^. 

We shall employ a now standard method of matched asymptotic expansions [52, 54-56, 65] to 
approximate u*' as e ^ 0, and use (18) and (22) to compute the asymptotics of ^^5 and A^^. Our 
treatment is facilitated by the fact that \f^ -^ Aq'^ = exponentially rapidly; up to exponentially 
small relative errors, it suffices when approximating v^^ in various regions of Q. to view it as a 
solution of £*uj^ = 0, rather than of £*uj^ = Af ^u*'. The alternative approach of approximating Xf' 
and v1 simultaneously is also possible, and has been applied to several models [11,51]. 

The following three asymptotic regions may be distinguished. Recall that we are considering 
the case when W attains a global minimum on dQ. at a distinguished saddle point H, so the limiting 
exit location a;* of the last section equals H. 

• A local (stable) region of size 0{e^l^) centered on S. 

• A region including all of ^ except the local region near S, and the boundary layer of width 
0{e^l^) near the characteristic boundary dO.. This region contains the fluctuational paths 
from the vicinity of S to the vicinity of H. 

• A boundary region, of size as yet unspecified but centered on H, and lying within the 
boundary layer. 

We now explain the corresponding asymptotic approximations for v^^: a Gaussian approximation 
near S, an outer approximation in the body of ^, and an inner approximation in the boundary 
region near H. The last, which we treat in Sections 5, 6 and 7, is where our treatment differs most 
radically from previous work; it will determine the asymptotics of the exit location density p^ 5. 
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4 The Outer Approximation 

The outer approximation to the quasistationary density v^ is a WKB expansion; equivalently, 
a characteristic (ray) expansion. We write 

v'^^(x) ~ K(x) exp {-W(x)/e) , e ^ 0, (23) 

for certain functions W : Cl ^ R and K : Q. ^ R normalized so that W(S) = and K(S) = 1, 
whose smoothness properties are as yet unspecified. We could instead attempt an approximation 
of higher order [68], with K(x) in (23) replaced by Ko(x) + eKi(x) + • • • . However we shall see 
in Section 9 that in the case of a characteristic boundary there may be difficulties with matching 
such outer expansions (in integer powers of e) to the inner approximation in the boundary region. 

Substituting the approximation (23) into the approximate forward Kolmogorov equation £*u'' = 0, 
and collating the coefficients of powers of e, yields equations for W and K: 






H{x\d,W) = (24) 

——(x% cW) + ^jrj7r^(x)j—-(x\ cm) K. (25) 



Here H : Q. xR^ ^ Ris the Hamiltonian (or energy) function 

H(x\ p,) = ^D'^(x)p,p, + b\x)p,, (26) 

and the eikonal equation (24) is the corresponding zero-energy Hamilton- Jacobi equation. For 
consistency, we have also used this Hamiltonian in the transport equation (25). The presence of 
a Hamilton- Jacobi equation suggests a classical-mechanical interpretation. In classical mechanics 
the Hamiltonian (26) would determine the motion of a particle onQ.; Q. x R^ would be interpreted 
as phase space, and p = (pi), i = 1, 2, as the momentum of the particle. But this Hamiltonian is 
precisely the Wentzell-Freidlin Hamiltonian governing the large fluctuations of the process x^(t) 
away from S [27] . Recall that if 



1 

r 



L{x\ i') = -D,j(x)(i' - b\x))(i' - Vix)) (27) 



is the Lagrangian canonically conjugate to H{x\pi), with the covariant tensor field D^j defined by 
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DijD^^ = 8i^ , then the Wentzell-Freidlin classical action function Vl^ : Q ^ ]R is defined by 

W{x)= inf / L(q,q)dt. (28) 

T>0 Jo 

q:[0,T]^n 
q(0)=5, q{T)=x 

That is, in Wentzell-Freidlin theory W(x) is computed as an infimum over trajectories from 5* to a;, 
and is a classical action function in the sense of classical mechanics. The Wentzell-Freidlin W will 
necessarily satisfy the Hamilton-Jacobi equation (24), so we may identify it with the function W 
in the outer approximation (23). 

By the calculus of variations, if the infimum in (28) is achieved by a trajectory extending from 
5* to a; , all portions of the trajectory that lie in Q. (rather than dQ.) must consist of least- action (zero- 
energy) classical trajectories [19]. Classical trajectories x(-) are the trajectories in Q. determined 
by H (or L); they satisfy the Euler- Lagrange equation 

d_ fdL\ dl^_^ 

Equivalently, if a classical trajectory is viewed as a pair of functions («(•), p(-)), specifying position 
and momentum (i.e., location in Q x M^) as a function of time, it must satisfy Hamilton's equations. 
Hamilton's equations are of the form 



/dB_\ 

dt[pj [-1 )\^ 

\dpk) 



(30) 



and express the fact that t i-^ (x(t),p(t)) must be an integral curve of a vector field on Q x M^ 
determined by H. For the Hamiltonian (26) they reduce to 



i'' = ^— = D'\x)p, + h\x) (31) 

Opk 

Pk = -^"T = dkD'^(x)p,pj + dkb\x)p,. (32) 

Notice that at all points along the trajectory, momentum and velocity uniquely determine each 
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other; by (31), 

p, = D,,{x)W - Vix)l (33) 

This is familiar from classical mechanics. 

If the least-action trajectory from 5* to a; exists and is unique, it is interpreted in Wentzell- 
Freidlin theory as the most probable (in the e ^ limit) fluctuational path from 5* to a;: the 
one whose cost is least. Since p^ = d^W at all points along the trajectory, W{x) (the cost of the 
trajectory) can be computed from 

W{x)= r p^dx\ (34) 

Js 

the line integral being taken along the trajectory. Equivalently, W satisfies 

W = p,i' = pdD''(x)p, + b\x)], (35) 

an ordinary differential equation which may be integrated along the corresponding trajectory in 
phase space. Note however that the infimum in (28) will not actually be achieved at finite transit 
time. It is readily verified for the Hamiltonian (26) that zero-energy trajectories emanating from 
a fixed point of the drift field b, such as S, will have infinite transit time. The formally most 
probable fluctuational paths require an infinite amount of time to emerge from S, and are more 
naturally parametrized by t G (— oo,r]. We discuss the physical consequences of this phenomenon 
elsewhere [56]. The integration of (35), as well as that of (31) and (32), must begin at t = — oo. 

For each such path limt^_oo x{t) = 5*, andby (33),limt^_ooP(t) = as well. But by examina- 
tion (S^ 0) is a fixed point (of hyperbolic type) of the deterministic flow on the phase space Q x M^ 
defined by Hamilton's equations. So in the language of dynamical systems, the most probable 
fluctuational paths of Wentzell-Freidlin theory may be viewed as forming the unstable manifold of 
the point (S^O) G H x M^. This manifold, which we shall denote A^(50)' i^ 2-dimensional; it is 
coordinatized (at least in a neighborhood of (S^ 0)) by (i) the choice of outgoing trajectory, and 
(ii) the arc length along the trajectory. M.'^gQ^ is a Lagrangian manifold: it is invariant under the 
flow. The stable manifold M^sfi) of ('"^^ 0) in Q x M^ is clearly the manifold p = 0, which by (31) 
and (32) is also invariant under the flow. It comprises 'cost- free' (AW = 0) trajectories that satisfy 
i' = b\x), and simply follow the mean drift toward S. 

The analogy with Hamiltonian dynamics can be pushed much further. The action W as com- 
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puted by (34) is a single-valued function of position on M'^soy ^^^^ Pi equals diW(x) at every 
point (x,p) G M'^sQ)- But it does not follow that W may be viewed as a single- valued function 
on Q.. This is because the unstable manifold may fold back on itself, and the projection (x^p) \-^ x 
from A^(5 0) to O. may not be one-to-one. Equivalently the zero-energy classical trajectories em- 
anating from S may cross each other, creating caustics [12, 14, 24, 36, 54-56]; these caustics are 
the projections onto Ci. of the folds of M'^soy ^^ shall not pursue the phenomenon of caustics 
at any length; unless otherwise stated we assume that the solution of the Hamilton- Jacobi equation 
is single-valued on Cl. In this case the map x \-^ p(x) will be single-valued. The projection 
map from M.'^^sq) to ^ may also fail to be onto. If this occurs, W(x), for at least some x ^ Q., 
cannot be viewed as the action of a classical trajectory from 5* to a; which lies entirely in Q.; the 
action-minimizing trajectory, if it exists, must pass through dO. on its way from 5* to a;. We defer 
discussion of this possibility to the next section. 

If we wish to compute the value of W at some point in Q. other than S, we can integrate 
Hamilton's equations from (S^O), continually updating x and p, until the specified point is reached. 
The ordinary differential equation (35), which may be integrated simultaneously, will yield the 
value of W at the endpoint. However the numerical computation of K, which has no elementary 
classical-mechanical interpretation, is more subtle. Since x' = dHjdpi, (dH/dpi)diK equals K, 
the time derivative of K along the trajectory emerging from S. So (25) becomes 

This equation cannot be integrated numerically unless the Hessian matrix didjW is known at all 
points along the trajectory. Fortunately d^d^W itself satisfies an ordinary differential equation, 
obtained as follows. Differentiating the Hamilton- Jacobi equation H = twice with respect to 
position on A^"5 0) yields 

/_d_ ^ dp^_d_\ /_d_ ^ dp^_d_\ ^ ^ ^ ^^^^ 

\dx^ dxi dpk / \dxi dx^ dpi j 
By rearranging terms, and using Hamilton's equations and d^p^ = d^d^W, we obtain 

(For notational convenience W^i^ signifies d^d^W henceforth; we also write h\j for djh\ etc.) 
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This version of the inhomogeneous Riccati equation for the Hessian W^ij is due to the present 
authors [55], though Ludwig apparently made use of an equivalent equation in the 1970's [49]. 
In any event, equation (38) may be integrated to yield W^ij at all points along an outgoing Wentzell- 
Freidlin trajectory. In numerical work the system (31), (32), (35), (36), (38) of coupled ordinary 
differential equations would be integrated simultaneously to produce W and K. 

At any x ^ Q., W^ij(x) = djPiix) specifies a 2-dimensional subspace of M"^, the tangent space 
T{x^p{x))M'l^sfi) to the unstable manifold ^^"50) C H x 1^ at («,/>(«)). The Riccati equation (38) 
can be viewed as determining a trajectory through the (compact) Grassman manifold of such 
2-dimensional subspaces. This sort of interpretation is standard in the theory of matrix Riccati 
equations [71], and facilitates the integration of (38) through points where W^ij diverges. Such 
divergences occur however only when the classical trajectory encounters a caustic [24]. This is 
because W^i^ diverges only at points x where the tangent space T(x,p{x))Mls o) 'turns vertical' [22]. 
We shall not pursue the consequences of caustics further here. 

Much more could be said about the geometric interpretation of the above system of equations, 
which is ultimately made possible by the symplectic structure of classical mechanics on Q x 1?. 
We confine ourselves here to pointing out the differential-geometric (coordinate-free) interpretation 
of the Wentzell-Freidlin fluctuational paths on Q.. The contravariant tensor field D'^ is naturally 
viewed as a Riemannian metric on ^, and may be used to raise and lower indices. If connection 
coefficients (Christoffel symbols) V jk are defined by 

r-. = iD-(D„,.A,,-ft„,) (39) 

in the usual way, the covariant derivative u'y of a vector field u' on Q. will be given by 

u' ,, = u\, +r ,ku^ . (40) 

It is easy to check that if the mean drift 6 = 0, the Euler- Lagrange equations (29) for the velocity 
field x\-) of the fluctuational paths on Q. reduce to the single covariant equation i'yi^ = 0. This is 
a statement of covariant constancy of the velocity of each such path along itself: the most probable 
fluctuational paths in this case are simply the geodesies of D^^ which emanate from the point S. 
If 6 7^ the situation is more complicated; computation yields 

{i' - b'hi' + d' - V)hf' = 0. (41) 
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In a similar way, if 6 = the Riccati equation (38), if viewed as applying to a function W 
defined on Q rather than on M.'^gQ^ simplifies to yield a covariant equation closely related to the 
equations of geodesic deviation on Q.. Taking 6 7^ yields a generalization. A fuller discussion of 
the differential-geometric interpretation, which owes much to work of Graham [31], may appear 
elsewhere. 

5 The Approximations Near the Fixed Points 

The outer approximation of the last section to the principal eigenfunction v*^ of C* must be 
supplemented by an inner approximation in a boundary region centered on H, and a Gaussian 
approximation in the stable region of size (9(e^/^) near S. To a certain extent it is possible to treat 
the approximations near H and S in parallel. 

This is possible because (H, 0), as well as (S, 0), is a fixed point (of hyperbolic type) of the 
deterministic flow on the phase space Q x M^ specified by Hamilton's equations (30), and the 
Wentzell-Freidlin Hamiltonian (26). In fact to any fixed point of the drift field 6 on Q (or Q.) there 
corresponds a fixed point of the flow on ^ x M^ (or Ci. x M^), at zero momentum. By examination, 
zero-energy classical trajectories which are incident on such points in phase space, as well as 
those which emanate from them, will have infinite transit time. As a consequence the zero-energy 
trajectories incident on (H, 0) may be viewed as forming the stable manifold of (H, 0) in Q x E^, 
which we shall denote M^fjQy 

We are assuming that the classical action W, computed from the variational definition (28), 
attains a minimum on dQ. at H. We shall also assume that the infimum in the expression (28) 
for W(H) is achieved; in particular, that W(H) is the action of a unique zero-energy classical 
trajectory q* : K ^ Q. which emanates from S at time t = —00 and is incident on i7 at t = 00. 
This trajectory q*, which is the formally most probable fluctuational path from S to dQ. as e ^ 0, is 
called the most probable exit path (MPEP). It may be viewed as a trajectory in phase space, in which 
case it necessarily lies in both the unstable manifold M'^sfi) ^'^d the stable manifold MljjQy The 
intersection M.'^gQ^ fl A^(/jo) of these two 2-dimensional manifolds will consist of zero-energy 
trajectories from (S*, 0) to (i7, 0). Examples are known in which the intersection consists of more 
than a single trajectory [55], but generically we expect that there is a unique zero-energy trajectory 
from S to H with minimum action. 

The behavior of W near S and H will constrain the possible approximations to u*' near S and H. 
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Since at all points on M'/^so) momentum p^ equals diW and the MPEP approaches («,/>) = (S^ 0) 
as t ^ — oo and (i7, 0) as t ^ oo, we expect that S and H are critical points of W on Cl. So to 
leading order, a quadratic approximation to W would seemingly be appropriate near both S and H. 
We shall shortly see that generically, the behavior of W near H (though not near S) is more 
complicated. But an exploration of the possible quadratic approximations to W will prove useful 
nonetheless. 

The assumption of quadratic behavior near the fixed points has the following consequences. 
Ast^— ooort^oo the left-hand side of the matrix Riccati equation (38) will tend to zero, 
yielding the algebraic Riccati equation 

^'"^ W,.,W. + ^^,. + ^W„, + ^ = (42) 



dpkdpi ' ' dx^dpk ' dx^dpk ' dx^dx^ 

for the Hessian matrix W^ij = W^ij(p) of second derivatives of W at the fixed point p, p = 5*, H. 
Here the partial derivatives oftheHamiltonian must be evaluated at («,/>) = (p,0). Substituting the 
explicit form (26) of the Wentzell-Freidlin Hamiltonian -ff (• , •) into (42) yields the matrix equation 

W,,D^'W,ki + W,,B^i + B^,W,,i = (43) 

for W^ij(p). Here we have written B = (B'j) for the linearized drift field b\j(p) at p, and D'^ 
signifies the local diffusivity tensor D'^(p). Recall that by assumption the matrix D'^ is a positive 
definite quadratic form at both S and H. We are also assuming that the matrix B'j is nonsingular 
at both S and H; this is the case for any generic drift field b. Note that equation (43) could be 
derived directly from the forward Kolmogorov equation [49, 57], without using the matrix Riccati 
equation (38). 

In matrix form equation (43) reads 

ZDZ + ZB + B'Z = (44) 

where Z = Z(p) = (W^ij(p)) is the 2-by-2 Hessian matrix of 1^ at p. Since W is computed 
from zero-energy trajectories emanating from S at t = — oo, the matrix W^ij(p) = dpi/dx^(p) 
specifies a 2-dimensional subspace of M"^, the tangent space T(^,o)M.'^g q-j of the unstable manifold 
^(3,0) ^t (p, 0). By assumption A^"5 0) extends from S to H along the MPEP q*,sop = H as well 
as p = S* is meaningful here. 
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We now digress to discuss the solution of the algebraic Riccati equations for the 2-by-2 
matrices Z(S) and Z(H). The analysis of the Riccati equation (44) is facilitated by a one-to-one 
correspondence between solutions Z (in general complex) and certain linear subspaces of C"^. 
This correspondence can be explained without reference to classical dynamics; it was in fact first 
explored in the context of optimal control theory [44, 70] . At each of the two fixed points p, define 
a 4-by-4 matrix T = T(p) by 




(45) 



T may be written as 



/ d^H 



d^H \ 



dpidx^ 



dp^dpj 



\ dx^dx^ dx^dpj / 




dx^dx^ dx^dpj 
\ dpidx^ dpidpj / 



(46) 



to reveal its symplectic structure; the representation (46) will be useful later. It is easily checked 
that a 2-by-2 (complex) matrix Z is a solution of the Riccati equation (44) if and only if the column 
space (over C) of the 4-by-2 matrix 



R 




(47) 



which is a 2-dimensional subspace of C^, is T-invariant, i.e., if and only if TR = RL for some 
2-by-2 complex matrix L. A 2-dimensional subspace of C^ which can be viewed as the column 
space of such a 4-by-2 matrix R (equivalently, a 2-dimensional subspace which is complementary 
to the subspace spanned by (0, 0, 1, 0) and (0, 0, 0, 1)) is known as a graph subspace. The one-to- 
one correspondence between solutions and subspaces facilitates the solving of the algebraic Riccati 
equation: one need only enumerate the 2-dimensional T-invariant subspaces of C"^, and determine 
which of them are graph subspaces. Each such graph subspace yields a solution Z, and these are 
the only solutions. 

The situation here is slightly more complicated in that for Z(p) we are interested only in 
solution matrices Z which are real symmetric. Gohberg, Lancaster and Rodman [30, 44] and 
Shayman [70] show in arbitrary dimensionality ('2' and '4' replaced by n and 2n respectively. 
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and T a 2n-by-2n real matrix) that if D > every n-dimensional T-invariant subspace of C^" 
is a graph subspace, and that Hermitian solutions Z of (44) will correspond to n-dimensional 
T-invariant subspaces M C C^" which are J -neutral in the sense that 



v^Jv = (48) 



for all vectors v ^ J\f . Here 




(49) 



is the fundamental symplectic form. The problem of finding real symmetric solutions Z thus 
reduces to constructing, from the spectral subspaces of T, n-dimensional J-neutral T-invariant 
subspaces whose corresponding matrices Z are real. 

This construction is much easier than it sounds. Since T is real, its eigenvalue set will be 
symmetric about the real axis, and since its eigenvalues are those of B and — S* the eigenvalue set 
will be symmetric about the imaginary axis as well. If all eigenvalues of T have unit multiplicity 
Lancaster and Rodman [44] show that there are 2^ possible subspaces M satisfying the conditions, 
where k is the cardinality of A: the set of eigenvalues A of T satisfying Re A > 0, Im A > 0. 
If no eigenvalue of T is pure imaginary, the possible M are formed from the 2^ subsets A' C A 
as follows: each is the spectral subspace of T corresponding to 

1. the eigenvalues in A', 

2. their complex conjugates, if distinct, 

3. the negatives of the eigenvalues in A \ A', and 

4. their complex conjugates, if distinct. 

If T has any pure imaginary eigenvalues, the construction of the 2^ possible M will be slightly 
different; each M must be expanded to include a unique subspace constructed from their eigenvec- 
tors [70] . If any eigenvalue of T has nontrivial multiplicity the set of possible subspaces M will no 
longer have cardinality 2^ . Nontrivial geometric multiplicity will engender an infinite number of 
subspaces J\f (and real symmetric solutions Z), though nontrivial algebraic multiplicity will yield 
only a finite (increased) number of real symmetric solutions [30, 44]. 
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Let us consider the implications of these results for the algebraic Riccati equations arising in 
the stochastic exit problem, which have n = 2. The eigenvalues of T are those of B, together 
with their negatives. But we are assuming that at both p = 5* and p = H, B is nonsingular. 
Since S is stable, either B(S) has two distinct eigenvalues Xs[(S), XsiiS) with negative real parts 
(both pure real, or Xsi(S) = XsiiS)*), or a single negative eigenvalue with nontrivial multiplicity 
(either algebraic or geometric). And since the smooth boundary dO. is characteristic, B(H) must 
have one negative eigenvalue Xs(H) < (whose eigenvector points along dQ.) and one positive 
eigenvalue A„(i7) > (whose eigenvector points into CI). At neither fixed point p does B(p) 
or r(p) have any pure imaginary eigenvalues. 

We see, accordingly, that at p = H there are normally k = 2 distinct (real) eigenvalues 
of r(p) satisfying Re A > 0, Im A > 0; namely, —Xs(H) and Xu(H). This will be the case unless 
— Xs(H) = XuiH); equivalently, unless trB(H) = b\i(H) = 0. We now explicitly rule out this 
possibility: we assume that —Xs(H) ^ A„(i7). This is not a major restriction, since drift fields b 
satisfying dib'(H) = are nongeneric. Recall that we have already defined the parameter 

fi = \XAH)\/X^(H). (50) 

The nongeneric case // = 1 will be seen below to be a difficult 'boundary' case, intermediate 
between the radically different cases /i < land// > 1 (the subjects ofSections 7 and 8 respectively). 
The difficulty with // = 1 is closely connected to the existence, by the results of Gohberg, Lancaster, 
Rodman, and Shayman, of an infinite number of possible solutions Z(H). This indeterminacy 
does not exist when // 7^ 1 ; there are precisely l'' = 4 real symmetric solutions of the algebraic 
Riccati equation (44) at p = H. We shall examine these four solutions shortly. 

We first dispose of the case p = 5*, which is comparatively straightforward. Here k = 2 except 
when B(S) has only a single (negative) eigenvalue, so one expects four solutions for Z(S) except 
in cases of nontrivial multiplicity. However it is easy to see that irrespective of the multiplicity, 
there will be only a single solution Z(S) of full rank. This is because if Z is of full rank, Z~^ exists, 
and multiplying (44) left and right by Z~^ yields 

D + BZ-^ + Z-^B' = 0. (51) 

This is a system of three linear equations for the three independent elements of Z~^. A bit of linear 
algebra shows that in the 2-by-2 case, the stability of B = B(S) will guarantee the existence of a 
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unique solution Z~^ and hence a unique nonsingular Z{S^. 

At p = S* the solutions to (44) of less than full rank, which include Z = 0, can be ruled out on 
physical grounds; S must be a local minimum of W, and the cost W(x) of fluctuations from 5* to a; 
must increase quadratically as a; — 5* increases. So the correct quadratic behavior of W near S is 
specified by the nonsingular solution. The nonsingular Hessian matrix W^ij(S) immediately yields 
an approximation to the eigenfunction u*' in the stable region of size (9(e^/^) near S. Since 

v'^^ix) ~ K(x) exp {-W(x)/e) (52) 

is the outer approximation in the body of ^, the natural approximation to v'^ near S is 

v'^^ix) ~ exp [-(x' - S'W,^J(S)(x' - S')/2e] . (53) 

This bivariate Gaussian approximation will match to the outer approximation if K is normalized 
so that KiS) = 1. Equivalently we could normalize the Gaussian so that it has total mass 
unity, and choose K(S) accordingly. In any event the covariance matrix Z(S)~^ of the Gaussian, 
which governs the small fluctuations of the process x^(t) on the (9(e^/^) lengthscale near S, 
can be computed numerically from (51); this has been known since the work of Ludwig [49]. 
On the (9(e^/^) lengthscale there are no well-defined most probable fluctuational paths: the process 
resembles a two-dimensional Ornstein-Uhlenbeck process. 

We now return to an analysis of the quadratic (or putatively quadratic) behavior of W near the 
saddle point H. A correct analysis depends crucially on the classical-mechanical interpretation of 
the matrix T. It follows from the representation (46) that linearizing Hamilton's equations (30) 
near (x,p) = (p, 0) yields 



(54) 



where ((5 a;, (5/>) = (a;,/>) — (p,0). In other words, r(p) specifies the linearization of the Hamiltonian 
flow near the fixed point (p, 0) in phase space. 

So the eigenvalues and eigenvectors of T(H) will have an interpretation in terms of the flow 
of zero-energy classical trajectories in the vicinity of H. We introduce some notation: let the 
eigenvectors of the linearized drift B(H) with eigenvalues Xs(H) < and Xu(H) > be denoted 
Bs and e„. By the definition (45), the eigenvectors of T(H) with eigenvalues Xs(H) and Xu(H) 



25 





will be (e^, 0) and (e„, 0). The eigenvectors of T(H) with eigenvalues —Xs(H) and — A„(i7) will 
be denoted (e„,^„) and (6^,^^); the interchange of subscripts is justified because —Xs(H) > 
and —Xu(H) < 0, indicating instability and stability respectively. Since (es,0) and (eg^gj are 
stable directions in phase space, the tangent space r(//o)A^(//o) of the stable manifold M^ffQ-j 
at (H, 0) will be their linear span (over E); similarly T(Hfi)M'(jj o) will be the linear span of (e„, 0) 
and(e„,^„). 

If W is determined in the vicinity of H by the zero-energy classical trajectories emanating from S 
and extending along M'^spy W,i]iH) = djPi(H) may be interpreted as specifying the tangent space 
T(Hfi)M.'(^SQ) of ^(5 0) ^t (-ff, 0). The four possible choices for the Hessian W^ij(H) found by 
solving the algebraic Riccati equation correspond to four possible choices for this 2-dimensional 
subspace of M"^. Each such corresponding subspace is spanned by the vectors (1,0, l^,ii, W^2i) 
and (0, 1 , 1^,12, 1^,22)- In other words T^Hfi)M}'s o) i^ the column space (over E) of the 4-by-2 matrix 



R(H) ={ . (55) 




This is precisely the restriction to M"^ of the subspace J\f used in the construction of Z(H). We found 
four (and only four) r(i7)-invariant J-neutral 2-dimensional subspaces J\f C C^ which yield real 
symmetric matrices Z(H), and we now see that their restrictions to IR"^ may be interpreted as the 
four possible choices for the tangent space T(^h,o)M'^s,o)- 

The interpretation of the four subspaces A/" as spectral subspaces of T(H) allows us to interpret 
the four possibilities for T(^h,o)M'(so) ^^ terms of classical dynamics. There are k = 2 positive 
eigenvalues of T(H), Xu(H) and —Xs(H). By the above construction, the four subspaces J\f will 
be the four spectral subspaces corresponding to the two eigenvalues ±A„(i7) and ±Xs(H). The 
corresponding possibilities for T(H,o)M.'^g q) are immediately expressible in terms of the eigenvectors 
of T(H). In particular the choice (-1-, — ) for the ± signs multiplying the two eigenvalues Xu(H) 
and Xs(H) yields T(^h,o)M.^hqj, and (—,-1-) yields r(//o)A^(//,o)- The two remaining possibilities 
are (-1-, -1-), which yields the linear span of (e„, 0) and (e^, 0), and (—,—), which yields the linear span 
of (e„,^„) and(es^gj. It is readily verified that the (—, — ) subspace corresponds to a matrix Z(i7) 
of full rank; the other three subspaces correspond to possible choices for Z(H) which are of less 
than full rank. The (-1-, -1-) subspace corresponds to the solution Z(H) = 0. 

For any stochastic model of the form that we are considering, the t ^ 00 limit of W^ij along 
the MPEP q* should exist, and equal one of the four possibilities for W^ijiH); equivalently, the 
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tangent space T(^x,p)M.'(s q) should converge to one of the corresponding possibilities for T(^h,o)M.'^s q) 
as («,/>) -^ (H,0) along the MPEP. The (+, — ) possibility, i.e., T^Hfl)M'^fjQy can be ruled out 
immediately. Since q* lies in A^"5 0) ^^^ by assumption q*(t) -^ H as t ^ oo, T^h,o)M}'sq-j must 
include a stable direction. But Xu(H) and —Xs(H) are both positive, so T(H,o)M]^ffQ^ cannot serve. 
The (— , +) possibility, i.e., T(H,a)Mljj q^ can also be ruled out. Since — A„(i7) and \s{H) are both 
negative, if T'(//o)A^(5 0) equalled r(//o)A^(//o) there would be an infinite number of zero-energy 
classical trajectories extending from (S*, 0) to (i7, 0), all with the same action. Such a situation 
would be highly nongeneric. 

There remain the (+, +) and (— , — ) possibilities for T{H,o)MlsQy ■^{+,+) ^^d M{--), as we shall 
call them, are 2-dimensional r(i7)-invariant subspaces of WJ^, spectral subspaces corresponding 
respectively to Xu(H), Xs(H), and to —Xu(H), —Xs(H). Since the stable eigenvalues of T(H) 
restricted to A/'(+,+) and A/'(_,_) are Xs(H) and — A„(i7) respectively, the MPEP q* will approach H 
exponentially, as exp(-\X,s(H)\t) if T(h,0)MIs^q) = J^{+,+) and as exp(-A„(i7)t) if T(h,o)MIs^o) = 
A/'(_,_). In principle either A/'(+,+) or A/'(_,_) could serve as T^H,0)M'^sQy but the relative magnitude 
of \Xs(H)\ and Xu(H) (i.e., the parameter //, or the sign of b\i(H)) turns out to determine which is 
generic. 

To see why this is the case, note that the MPEP, viewed as a trajectory in A^"5Q) C H x M^, 
will approach (i7, 0) in a manner specified by the linearized Hamilton's equations (54), and that 
the negative eigenvalues of the linearized drift matrix T(H) in (54) are Xs(H) and —Xu(H). 
Generically, q* will approach (i7, 0) as t ^ oo along the less contractive direction in phase space; 
the stochastic model would have to be carefully 'tuned' to arrange for the incoming MPEP to 
approach along the more contractive direction in (54). If // < 1, —Xu(H) < Xs(H) < and 
Xs(H) is less contractive; in this case the MPEP will generically approach H as exp(—\Xs(H)\t), so 
if /i < 1 then T^Hfl)M]^sQ) = J^{+,+) will be generic. Similarly if // > 1, Xs(H) < — A„(i7) < and 
the MPEP will generically approach H as exp(— A„(i7)t), so if // > 1 then T^Hfl)M]^sQ^ = A/'(_,_) 
will be generic. Consequently, in phase space the MPEP will generically approach the fixed 
point (i7, 0) along the tangent vector (e^, 0) if b\i(H) > 0, and along (eg^gj if b\i(H) < 0. 

We now consider the extent to which the behavior of W near H can be quadratic. A great deal 
of light is thrown on this question by a sketch of the pattern of zero-energy classical trajectories 
emanating from S, when prolonged to the vicinity of H. We may regard these trajectories, which 
include the MPEP g* , as lying in the manifold M "50)- But in a sufficiently small neighborhood of H, 
the flow of these trajectories on Mfs o) should be given to high accuracy by the the corresponding 
flow on the tangent space T(H,o)M}'sQy If fi < 1, this tangent space is generically A/'(+,+), the linear 
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Figure 2: The generic appearance of the flow near H of the zero-energy classical trajectories 
emanating from S; equivalently, the Hamiltonian flow on the manifold M.'^gQ-^ projected 'down' to 
configuration space by the map (x,p) \-^ x. Part (a) of the figure illustrates the case ji < 1, i.e., 
dib\H) > 0; part (b), the case /i > l,i.e., dib\H) < 0. In both cases the MPEP [most probable exit 
path] is the solid curve incident on H. The wedge-shaped shaded regions are classically forbidden; 
also, the dashed trajectories (which extend beyond the e^ ray, which lies in dQ.) are unphysical. 

span of the stable vector (e^, 0) and the unstable vector (e„, 0); if // > 1, it is generically A/'(_,_), 
the linear span of the stable vector (Ss^gJ and the unstable vector (eu^gj. It follows that the 
zero-energy classical trajectories emanating from S, if projected 'down' from M'^so) to ^ by the 
map («,/>) 1-^ X, will generically trace out in the immediate vicinity of H the hyperbolic patterns 
shown in Figures 2(a) and 2(b). The principal axes of the hyperbolae are e^ and e„ if // < 1, and 
Eg and e„ if /i > 1. By convention, we assume that the e^ ray is vertical, and that the MPEP 
approaches H from the first quadrant. 

Figure 2(a) reveals a phenomenon not previously suspected: since the stable eigenvector e^ 
of B(H) lies in d^, the MPEP q* will generically be tangent to dD. for stochastic models with // < 1 . 
This prediction of a 'grazing MPEP' has been confirmed by numerical studies. Also, since the 
(-I-, -I-) solution for Z(H) is the 2-by-2 zero matrix, when // < 1 one expects that that W^ij(x) -^ 
as X ^ H along the MPEP. This too has been confirmed by numerical studies. The consequent 
'flat' behavior of W near the approaching MPEP has an appealing physical explanation: near H 
the most probable fluctuational paths (perturbations of the grazing MPEP) come near to following 
the drift field 6, so very little additional action (cost) is built up as H is approached. 

We have not commented yet on the most visible feature of Figures 2(a) and 2(b), which is 
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highly disconcerting: for both // < 1 and // > 1 there is a wedge-shaped region near H, beyond 
the boundary ray e„ and the boundary ray e„ respectively, which is not reached by any of the 
fluctuational paths. This wedge, the generic presence of which has been confirmed numerically, is 
'classically forbidden' in that it cannot be reached by any zero-energy classical trajectory emanating 
from S. Unless the boundary ray lies in dO., the wedge has nonempty interior. The generic 
presence of an unreachable region has a very simple interpretation: Generically, the projection 
map (x,p) 1-^ X from M^sfl) to Q. is not onto; its range does not include all ofQ.. The possibility 
that the projection map might fail to be onto on account of folds in the manifold A^(5 0) ^^^ 
mentioned in Section 4, but that this failure is generic is a new result. 

We expect that Wix), for any x in the wedge, is the action of an action-minimizing trajectory 
which passes through dO. on its way from 5* to a; [19]. Only the portions of the trajectory that 
lie in Q. (rather than dO) will be classical, i.e., will be solutions of the Euler- Lagrange equation. 
The trajectory will be piecewise classical: it will have at least one 'corner,' or bend. The fact that 
generically, such fluctuational paths need to be considered for at least some endpoints x in any 
neighborhood of H has not previously been realized. 

By a careful optimization over piecewise classical trajectories it is possible to work out, for 
X in the wedge, the leading dependence of W(x) on 6x = x — H in the vicinity of H. But we 
may economize on effort by applying our previous results. Suppose for the moment that to leading 
order W(x) in the wedge behaves quadratically as x ^ H. Then the quadratic dependence will be 
specified by a limiting Hessian matrix Z(H) = (W^ij(H)), which must satisfy the algebraic Riccati 
equation (44). There are exactly four possibilities for Z{H)\ equivalently, four possibilities for 
the corresponding graph subspace J\(H) C R^. But T(Hfi)M^fj,Q)' ^•^•' sp{(e„, 0), (e„,^„)}, is the 
only one of the four possible subspaces on which the linearized Hamiltonian flow (45) has positive 
eigenvalues: Xu(H) and —Xs(H). The other three have at least one stable direction, and would 
correspond to a pattern of action-minimizing trajectories in the wedge which converge on H, rather 
than spreading outward from H. This is impossible, so we must have J\(H) = T^Hfl)M'^fj q)- If ^ 
in the wedge behaves quadratically near H then the limiting Hessian matrix Z(H) will necessarily 
be the rank-1 2-by-2 matrix corresponding to T^H,0)M]^ffQy and will differ from Z(H). 

The extent to which the action-minimizing trajectories in ^ x M^ giving rise to this Z(H) are 
only piecewise classical follows from the fact that they must emanate from (i7, 0), and thereafter 
(near H) lie in the subspace of IR"^ spanned by (e„, 0) and (eu^gu), the two expanding eigenvectors 
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Figure 3: The generic appearance of the flow near H of the most probable fluctuational paths 
{i.e., the action-minimizing trajectories) emanating from S. The e^ ray lies in 9Q. Parts (a) and (b) 
of the figure illustrate Subcases A and B of the case // < 1; they differ in the relative placement 
of the e„ and e„ rays. Part (c) illustrates the case // > 1. In all three parts the MPEP [most 
probable exit path] is the solid curve incident on H. The other trajectories include zero-energy 
classical trajectories of the sort shown in Figure 2, and 'bent' trajectories which are prolongations 
(piecewise classical rather than classical) of the MPEP The classically forbidden wedge-shaped 
shaded regions are reached (as e ^ 0) only via bent trajectories. In part (a) these trajectories have 
only one bend, at H. But in parts (b) and (c) the dashed portions of the trajectories emanating 
from H are unphysical, and the trajectories entering the shaded regions have two bends: they 
extend along d^ before reentering Q.. 
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of T(H). Near H each such trajectory, when parametrized by t > — oo, should be of the form 

t^(i7,0) + C„exp(A„(i7)t)(e„,0) + C'„exp(|A,(i7)|t)(e„,^J, (56) 

for C'u, C'u two constants that determine the trajectory. In general these action-minimizing tra- 
jectories will be tangent (as t -^ — oo, or as the trajectories emerge from H) to whichever is the 
less expansive of the two expanding eigenvectors of T(H). If // < 1 this is (e„,^„); if // > 1 it 
is (e„, 0). This tangency condition determines the behavior near H of the velocity field x(-) of the 
action-minimizing trajectories in the wedge. If // < 1 the portions of the trajectories emanating 
from H should be tangent to e„. If // > 1 they should be tangent to e„. 

Figure 3 is an extended version of Figure 2, which includes the rays e„ and e„, and the 
additional trajectories (emanating from H, and tangent to them) predicted by the assumption of 
quadratic behavior of W in the wedge. The // < 1 case has two subcases, shown in Figures 
3(a) and 3(b). We shall refer to them as Subcase A and Subcase B; they differ (as regards the 
projected eigenvectors of T(H)) in only one essential way: in the relative placement of the rays 
e„ and e„. We shall show in the next section that when // > 1, the e„ ray necessarily lies between 
the Bg and e„ rays. This situation is shown in Figure 3(c); when // > 1, there are no subcases. 

Subcase A of the case // < 1 is the most straightforward. It is clear from a glance at Figure 3(a) 
that in this subcase, the action-minimizing trajectory from S to any point in the wedge is a 
prolongation of the MPEP: it extends from S to H, experiences a discontinuous change in direction, 
and then enters the wedge. Interestingly, points on dQ. that are on the 'wedge side' of H are reached 
only via trajectories that make a second passage through Q.. 

We stress that these 'bent' trajectories have a probabilistic interpretation. By Wentzell-Freidlin 
theory, any point x in the interior of the wedge is (as e ^ 0) preferentially reached during large 
fluctuations away from S* by a bent trajectory. That is, when // < 1 the fluctuation will in Subcase A 
preferentially drive the random proces s a; ^ (t ) to the vicinity of i7 via the MPEP q * , before the wedge 
is traversed and the vicinity of x is finally reached. Although Figure 3(a) displays only the portion 
of the wedge in the vicinity of H (i.e., the wedge as computed in the linear approximation), studies 
of particular models show that it may extend a considerable distance from H. The numerical 
computation of W(x), for x in the wedge, requires an integration along the appropriate bent 
trajectory terminating at x. 

There is a problem extending these conclusions to Subcase B, and to the case // > 1. It is clear 
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Case [Subcase] 


Z{H) 


Z{H) 


M{H) 


Ar(i7) 


fi< \,i.e.,dMiH)>Q [A] 





rank- 1 


A/'(+,+) = sp{(e„0),(e„,0)} 


T(Hfi)M'(Hfl) 


/i < 1, i.e., d,U{H) > [B] 







Ar(+,+) = sp{(e„0),(e„,0)} 




fi> I, i.e., d,U{H) < 


rank- 2 




A/'(-,-) = sp{(e„^,),(e„,^J} 





Table 1: The generic limiting behavior (as x ^ H) of the Hessian matrix Z(x) = {W^ij(x)), and 
of the corresponding graph subspace J\f(x) C M"^. ^(a;) ^ Z(H) if x ^ H from outside the 
wedge, i.e., along the manifold M'^soy Similarly A/'(a;) -^ J\f(H). In Subcase A of the case // < 1, 
Z(x) -^ Z(H) as X ^ H from within the wedge; similarly ^f(x) -^ Af{H). Each of Af{H) 
and Jv{H), if defined, is the linear span of a pair of eigenvectors of the linearized Hamiltonian 
flow T(H), since the tangent space T(h,o)M'^h,Q) equals sp{(e„,0), (e„,^„)}. In Subcase B and 
when /i > 1 the behavior of the action in the wedge is not quadratic near H. 

in Figures 3(b) and 3(c) that classical trajectories which are of the form 



t ^ i7-i-C„exp(A„(i7)t)e„-i-C'„exp(|A,(i7)|t)e, 



(57) 



and which are tangent to the e„ ray (resp. the e„ ray), are necessarily unphysical. They are 
the dashed trajectories, which penetrate into the complement of Q. before returning to dO. and 
passing through the wedge. This situation is quite different from Subcase A, where the trajectories 
emanating from H penetrate immediately into the wedge. We expect therefore that in Subcase B 
and when // > 1 the true action-minimizing trajectory from H to any point x in the wedge will 
extend along 9Q before reentering the wedge. Moreover, since the putative Hessian matrix Z{H) 
corresponds to an unphysical (impossible) velocity field x{-) for the most probable fluctuational 
paths, in Subcase B and when // < 1 we do not expect the action in the wedge to behave quadratically 
near H. The flow field of the action-minimizing trajectories within the wedge may differ from the 
predictions of the quadratic approximation. 

We summarize our conclusions in Table 1. The behavior of W near H, both inside and outside 
the wedge and for both // < 1 and ji > 1 , is quadratic and fully understood with the exception 
of the behavior inside the wedge in Subcase B, and when // > 1. Even these difficult cases lend 
themselves to a partial analysis. The precise behavior of W inside the wedge in Subcase B and 
when /i > 1, even in the linear approximation near H, can only be obtained by solving a difficult 
variational problem. But it is clear on physical grounds that W(x) approaches W(H) quadratically 
as X approaches H along dO. from the wedge side of H. This is because on the wedge side of H, 
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the cost function W(x) on dO. satisfies (in Subcase B, and when // > 1) 

W(x)-W(H)= r p^dx\ (58) 

Jh 

the line integral being taken along dO. from H to x. Here the momentum p = p(x), at any 
point X G do., takes a value determined by (i) the zero-energy constraint, and (ii) the constraint 
that the corresponding velocity x lie in dO.. 

In Subcase A, in general Z(H) ^ Z{H)\ in Subcase B and when // > 1 there is no reason 
why the quadratic growth of W along dQ., on the wedge side of H, should equal the quadratic 
growth on the other side (arising from Z(H), rather than from (58)). We conclude that generically, 
W fails to be twice continuously dijferentiable at H. As x ^ H from within ^, in Subcase A 
the Hessian matrix W^ij(H) will converge to Z(H) or Z(H) depending on whether x ^ H from 
outside or inside the wedge. The corresponding graph subspace J\f(x) will tend either to JV(H), 
i.e., to r(// 0)^(5 0)' o^ to JV(H). The situation for Subcase B and // > 1 is similar, except that 
we lack a precise description of the behavior of W inside the wedge. Notice that if s denotes arc 
length along dCi. measured from H, with s < on the wedge side, the value d^W/ds^(s = 0+) 
will equal e\Z{H)es, and the value d'^Wjds^is = 0—) will equal e\Z{H)es (in Subcase A), or a 
quantity derived from (58) (in Subcase B, or if // > 1). Generically these two one-sided second 
derivatives will be different; in fact if // < 1 the former, as we have noted in Table 1, will equal 
zero. This justifies the statements about d'^Wjds^is = 0±) made in the Introduction. 

It is worth noting that W(x) will generically fail to be a twice continuously differentiable 
function of x not merely at x = H, but along the entire boundary of the wedge. This is strongly 
reminiscent of the Stokes phenomenon, which occurs in the asymptotic expansion of the solution 
(in the complex plane) of an ordinary differential equation with a small parameter multiplying the 
term with highest derivative. The exponent of the dominant exponential factor in the expansion 
is non-smooth along certain curves ('anti-Stokes lines') emanating from turning points. The 
implications of the Stokes phenomenon for the asymptotic expansion of Airy functions, etc., are well 
known. It is probably best to view the phenomenon of the wedge as a two-dimensional counterpart 
to the conventional Stokes phenomenon. For the forward Kolmogorov equation on Q c M^, which 
is a partial differential equation, the anti-Stokes line appears as a curve in Q. (the boundary of the 
wedge) rather than in the complex plane. (Cf. Berry and Howls [5].) In order to see it, no analytic 
continuation is required. 

In the following three sections we approximate the principal eigenfunction v^^ and the pro- 
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cess x^(t) in the boundary region near H. Our results on the flow of classical (and piecewise 
classical) trajectories near H will to a large extent dictate the choice of asymptotic approximations. 

6 Matched Asymptotic Expansions, Continued 

It is clear, in broad outline, how to apply the method of matched asymptotic expansions sketched 
in Section 3. In Section 4 we derived an outer approximation to the principal eigenfunction uj^ 
which is valid in the body of Q.; it will necessarily match to the Gaussian approximation (53) in 
the stable region of size (9(e^/^) surrounding S. We lack only an inner approximation to u*' valid in 
the boundary region near H. If it is to match to the outer approximation, the inner approximation 
must in the far field have leading asymptotics 

v'^(x) ~ exp [-(x' - H')W,,j(H)(x' - H')/2€] . (59) 

But as we have seen from the last section, this expression must be interpreted with care. By con- 
vention W^ij(H) signifies the limiting Hessian matrix obtained by taking x ^ H along the most 
probable exit path (MPEP), or in general from within the complement of the classically forbidden 
'wedge.' The asymptotic statement (59) is interpreted as holding when (x — H)/e^^^ -^ oo within 
the complement of the wedge. If the behavior of W in the wedge is quadratic near H, the limiting 
Hessian matrix obtained by taking x ^ H from within the wedge is denoted W^ij(H). We have 
seen that in Subcase A of the case dih\H) > 0, the behavior of W in the wedge is indeed quadratic 
near H, so the counterpart to (59) (involving W^ij{H)) should hold as {x — H)/e^^^ -^ oo within 
the wedge. In other words the far-field asymptotics of the inner approximation to v^ must in this 
case he piecewise bivariate Gaussian: different decay (or growth) rates will be found in the wedge 
and its complement. In Subcase B or when dih\H) < (z.e., // > 1, in the notation of the last 
section), the behavior of W in the wedge will not necessarily be quadratic near H. But the inner 
approximation to u*' must still have Gaussian asymptotics as {x — H)/e^^^ -^ oo along d^, with 
different decay rates on the two sides of H. 

The inner approximation to u*' can be found, at least in principle, by solving a linearized version 
of the (approximate) forward Kolmogorov equation C*p = Oin the boundary region near H. In the 
linear approximation we take b'(x) sa B'j(H)(x^ — H^) and D'^(x) sa D'^(H), so the Kolmogorov 



34 



equation reduces to 

(e/2)D^^(H)d,d,p - d,[B^,(H)(x^ - W)p] = 0. (60) 

Assume for the moment that the appropriate lengthscale on which the inner approximation should 
be defined is indeed the (9(e^/^) lengthscale. If so, we can employ the 'stretched' variable X' = 
(x* — W)l e^l'^, in terms of which (60) becomes 



1 .. d^p d 



■^D'^H) ^^.'^^. - B' j(H)^^^[X' p] = 0. (61) 



We can change variables to reduce this covariant equation to a noncovariant, but more under- 
standable form. Under a linear change of variables (X'y = L'jX\ i.e., X' = LX, the matrices 
B(H) and D(H) transform to LB(H)L-^ and LD(H)L respectively. Choosing L = D{H)-^I^ 
transforms D(H) to the identity matrix. But since H is a saddle point, irrespective of coordinate 
transformations the linearized drift B(H) will have one positive eigenvalue (Xs(H)) and one neg- 
ative eigenvalue (A„(i7)). By a further change of variables (a rotation) we can arrange matters so 
that B^2(H) = 0, and 

B(H) = "' ' , (62) 

for some real constant c. The constant c is not determined by Xu(H) and Xs(H). Since D(H) = I 
is preserved under rotations, with respect to the new system of coordinates equation (61) becomes 

llBr-*\lBr- ^"(^•aJ^t^-Vl - UH>^.iry) - c^[A- VI = 0. (63) 

In terms of the transformed coordinates (X\X^) we may view the region Q. as the right-half 
plane X^ > 0, and its boundary dO. as the X^-axis. This was the convention of Figures 2 and 3. 

This system of coordinates is computationally easy to work with. Suppose for simplicity that 
Xu(H) = 1; this is an innocuous normalization condition that can be absorbed into a redefinition 
of time t (and noise strength e). Then 



B(H) = , (64) 
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where // = \Xs(H)\/Xu(H) as in the last section. The forward Kolmogorov equation (63) reduces 
to 

Moreover a bit of matrix computation, using the form (64) for B(H) and D(H) = I, yields 

(e„0) = (0,1,0,0) (66) 

(e„,0) = (/i+l,c,0,0) (67) 

(e,,^,) = (/i-l,c,-2/i + 2,0) (68) 

(^u,9u) = (-2/ic, /i^ - 1 - c^ -2/ic(/i - 1), 2/i(/i^ - 1)) (69) 

for the four eigenvectors of the linearized Hamiltonian flow T(H) at the point (i7, 0) in phase 
space, as discussed in the last section. (Normalization is irrelevant here; the negatives of these 
vectors could equally well have been chosen.) The fact that e^ = (1, 0) is in agreement with the 
convention of Figures 2 and 3. 

The formulae (66)-(69) explain the positioning of the rays e„, e^, and e„ in Figures 3(a), 
3(b), and 3(c). Recall that in those figures the MPEP is taken to approach from the first quadrant; 
if /i < 1, it is generically tangent to e^, and hence to the positive X^-axis. By examination 
of (66)-(69), if /i < 1 and c < then e„ will lie between the positive X^-axis and e„, while 
if c > then e„ (taken to point into the right half-plane) will lie between the positive X^-axis 
and e„. So the Subcases A and B of Figures 3(a) and 3(b) are simply the subcases c < and c > 0. 
This correspondence assumes of course that the MPEP is tangent to the positive X^-axis; if it 
approached H from the fourth quadrant, and were tangent to the negative X^-axis instead, then the 
interpretation in terms of sgn c would be reversed. 

Figure 3(c) is justified, as well. If // > 1 the approaching MPEP is known by the argument of 
the last section to be generically tangent to e^. By (68), e^ oc (// — 1, c), and our convention that the 
MPEP approaches from the first quadrant mandates that c > 0. It is easy to verify, by examining 
(66)-(69), that in the right-half plane, when // > 1 and c > the e„ ray necessarily lies between 
the the e^ ray and the the e„ ray. This justifies the relative positioning of the rays in Figure 3(c). 

Now that a canonical system of coordinates has been chosen, the rank-2 Hessian matrix 
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Z(H) = (W^tj(H)) may be computed; substituting (64) into (51) yields 

Recall that generically this limiting Hessian matrix arises (when a; ^ i7 in the complement of 
the wedge) only when // > 1. If // < 1, Z(x) -^ as x ^ H along the MPEP, and as noted in 
Table 1 we have Z(H) = instead. In Subcase A of // < 1, the limiting Hessian matrix Z(H) 
in the wedge exists, and equals by Table 1 the rank-1 matrix corresponding to the graph subspace 
J\f(H) = T(h^q)M'Ijjq^ in E"^. This subspace is the linear span of (e„,0) and (e„,^„), and the 
explicit formula 

for the associated solution of the algebraic Riccati equation follows by some elementary manip- 
ulations. In Subcase B, and when ji > 1, the behavior of W{x) in the wedge as a; ^ i7 is not 
expected to be quadratic. However it follows from the formula (58) that its behavior as x ^ H 
along d^ (from the wedge side) is quadratic, with limiting second derivative 

Za^(H) = 2fi. (72) 

This limiting second derivative has a simple interpretation: on the wedge side of H, the cost 
(action) of the most probable trajectory leading to any point on dO. arises from the drift b on dO. 
itself. 

The formulae for the one-sided second derivatives Z(H), Z(H), and Zgo^H) make quite precise 
the far-field {X -^ oo) Gaussian (or inverted Gaussian) asymptotics which must be imposed on the 
solution p(X\X^) of the transformed Kolmogorov equation (65). First, we must have to leading 
order 

p(X\X^) ~ exp [-X'Z(H)X/2) , X ^ oo outside the wedge. (73) 

Also, in Subcase A of // < 1 we must have 

p(X\X^) ~ exp [-X'Z(H)X/2) , X ^ oo inside the wedge, (74) 
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and in Subcase B of // < 1, or when ji > 1, we must have 



p{X\X') ~ exp -Z9a(i:^)(XY/2 



X = (0, X^) ^ oo, on the wedge side, (75) 



since the boundary dO. is the X^-axis. The location of the wedge here is as shown in Figure 2. 
If /i > 1 the wedge is the sector bounded by the rays e^ and e„ which contains neither e^ nor e„; 
if /i < 1 it is the sector bounded by the rays e^ and e„ which contains e„ but not e^. The Dirichlet 
(absorbing) boundary condition p(0, = must also be imposed. It expresses the fact that the 
quasistationary density is absorbed on the boundary. 

In general it is not easy to solve the partial differential equation (65) on the half-plane X^ > 0, 
subject to these boundary conditions. Our treatments of the generic // < 1 and generic // > leases, 
in Sections 7 and 8 respectively, are designed to circumvent this problem. For the case // < 1 
we shall expand on a larger lengthscale than (9(e^/^); for the case // > 1, on which we have less 
information, we shall use stochastic analysis. In advance of our detailed treatments, we observe 
that the far-field behavior of the exit location density p^ (on the 0(e^/^) lengthscale near H) follows 
from (73), (74), and (75). By (18) the exit location density is simply proportional to the normal 
derivative of the inner approximation uj^. If we substitute the known expressions for Z(H), Z(H), 
and ZdoiH), we obtain far-field asymptotics 



Pe{s) ~ < 



exp 

exp 



• (^Ve)] , 

pip + 1)^ 
c^ + {p + !)• 



is^le) 



sje'l^ 



s/e 



1/2 



-l-oo. 



-oo 



(76) 



if /i < 1 (Subcase A), and 



Pe(s) ~ < 



exp [O • (s7e)] , s/ei/2 ^ +oo 
exp — /i(s^/e) , s/e^/^ ^ — oc 



(77) 



if /i < 1 (Subcase B), and 



Pe(s) 



~ < 



exp 
exp 



pip - 1) 

~C^ + ip- 1)2 

-pi-s^/e) 



isVe) 






-l-oo, 



-oo 



(78) 



if p > 1 . The zero coefficients in the exponents of (76) and (77) indicate that the corresponding 
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asymptotics are sub-Gaussian. We have written s for Ax^ (the distance along dQ., measured 
from H), for consistency with the Introduction, and have taken the wedge side of H to be the side 
on which s < 0. Notice that as // | 1 (resp. /« i 1) the asymptotics of (77) and (78) come into 
agreement. Subcase A, however, has no // > 1 counterpart. 

There is one degenerate case in which the partial differential equation (65), subject to the 
boundary conditions (73)-(75), can be solved explicitly. This is the case when c = 0, i.e., when 
the linearized drift B(H) at the saddle point has the property that its eigenvalues e^ and e„ are 
orthogonal with respect to the inner product specified by the local diffusivity tensor D(H). If c = 
and /i > 1, the classically forbidden wedge vanishes: the rays e^ and e„, which form the boundary 
of the wedge, become identical. In picturesque language, as c ^ (when // > 1) the wedge 
disappears into the boundary dQ., and the action becomes twice continuously differentiable on a 
neighborhood of the saddle point. In this case it follows from (78), by setting c = 0, that the two 
Gaussian decay rates of the asymptotic exit location density are equal. This suggests that p^ is 
asymptotic to a Gaussian. The confirmation of this, however, requires the solution of (65). 

If there is no wedge, the imposed far-field asymptotics become Gaussian rather than piecewise 
Gaussian; they reduce to p(X^, X^) ~ exp (-X*Z(i7)X/2) , since only (73) applies. But if c = 
the expression (70) for the matrix Z(H), which is valid when // > 1, simplifies to 



Z(H) =11. (79) 




So the far-field asymptotics which must be imposed on the solution of (65) simplify greatly; we have 
p(X\X^) ~ exp[(XO^ — fi(X^)^]. Moreover if c = the final term in (65) vanishes, allowing a 
solution on X^ > with these asymptotics to be found by separation of variables. It is 

p(X\X^) = exp [-p(X^f] G(X') (80) 

in which 

G(z) = e^ erf(z) = ^e^' ^ e"^' ds. (81) 

a/tT Js=0 

When c = therefore (and p> 1) we have the inner approximation 

vl{x\x^) ~ {/ai7)exp (-W/(i7)/e)}exp \-p{x^ - E^fj^ G [{x^ - H^)/e^^^) (82) 
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in which we have included the prefactor K(H)exp {—W(H)/€) to facilitate matching with the 
outer approximation. 

Note that in (82) we must interpret x^ — H^,x^ — H^ as linearly transformed ('primed') versions 
of the original coordinates x^ — H\ x^ — H^. As remarked, X = (X\X^) in (80) really means 
X' = OD{H)~^I^X, for O a suitably chosen (orthogonal) rotation matrix; similarly, x — H 
in (82) must be interpreted as OD{H)~^l^{x — H). Using (79), we can rewrite (82) in a partially 
covariant form as 

v%x\ x^) ~ {A^(i7)exp (-W/(i7)/e)} exp [-{x' - W)W,,j{H){x' - H')/2e] erf [(x^ - H^)/e^^^] 

(83) 

with the same proviso on the interpretation of the factor x^ — H^ in the argument of the error 
function. This way of writing the inner approximation to u*' makes it clear that its leading 
asymptotics as x/e^^^ -^ oo are indeed those of (59). It also reveals that it is closely related to 
the Gaussian approximation (53), which is valid near S. The final factor in (83), which has no 
counterpart in (53), is attributable to the absorbing boundary condition at X^ = 0. This boundary 
layer function may be written in fully covariant form as 

erf (Xu(Hy/\'D,,(H)(x^ - H^)/e'/^) (84) 



where n' is a contravariant unit normal vector to dQ. at H, satisfying n'Dij(H)fJ = for any 
vector f tangent to dQ. at H, and normalized so that n'Dij(H)n^ = 1. We have included in the 
argument of the error function a factor XuiH)^^^, which will be present if Xu(H) ^ 1. 

We stress that we are able to find such a simple inner approximation as (83) in the case // > 1 
only when c = 0, i.e., only when the eigenvectors of B(H) are orthogonal with respect to the 
inner product specified by Dij(H). A covariant way of expressing this condition may be derived 
as follows. Bs'D.jeJ = 0, i.e., elDiHy^e^ = 0, means that D{H)-^l^e, and DiHy^'^Bu 
are orthogonal in the conventional sense. Equivalently, D(H)~^^^B(H)D(Hy^^ has orthogonal 
eigenvectors and is a symmetric matrix. But D(H)~^^^B(H)D(Hy^^ if and only if B(H)D(H) is 
symmetric, i.e., if the tensor h\jD^^ is symmetric under the interchange of i and katx = H. This 
is the covariant form of the condition. It is easily checked that if 6' = —D^Wj^^ for some scalar 
potential field $ which has a saddle point at H, then the condition is necessarily satisfied. The 
c = condition for the validity of the inner approximation (83), when // > 1, is really a condition 
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that the drift b be locally gradient at H. 

Stochastic exit models in which this local gradient condition holds are (sadly) nongeneric. But 
if /i > 1 , and the condition holds, it is an easy matter to apply the technique sketched in Section 3 to 
the inner approximation (83) to determine both the asymptotics of the exit location density on dO., 
and the asymptotics of the principal eigenvalue Xf^ of £*. Equation (18) yields 

p,(x) oc exp [-(x' - H')W,,j(H)(x^ - H')/2e] , x e dQ (85) 

for the density of the exit location distribution in the e ^ limit, on the 0(e^/^) lengthscale 
near H. So the density on dQ. is indeed asymptotically Gaussian, with the same falloff rate 
as (x — H)/e^^^ ^ oo to either side of the saddle point. Once again we stress that this Gaussian 
behavior can only occur in the absence of a classically forbidden wedge. 

Equation (22), which expresses Xf^ in terms of the flux of probability into dO. near H, when 
applied to the inner approximation (83) yields 



(Et,)-i ~ Af ~ -^det Z(S)J^;^K(H) exp {-W(H)/e) e ^ 0. (86) 

Here the factor A/det Z(S) arises from the denominator of (22). We have approximated the integral 
in the denominator of (22) by an integral of the Gaussian approximation (53) to v^ over the stable 



region of size (9(e^/^) near S; this equals IttJcJ det Z{S). The asymptotics of (86) are valid only 
if the coordinates x^ near H have been linearly transformed in such a way that D(H) = I; this 
was of course assumed in the derivation of the inner approximation (83). 



(Er,)-^ ~ Af ~ ^'^ iT.ffl Kimcxp {-W(H)/e) , . ^ (87) 

TT y I det Z/ {Jd ) I 

is the generalization to arbitrary coordinate systems. 

Wentzell-Freidlin theory provides the leading order growth of the MFPT Er^ as e ^ 0; it is 
exponential, with rate constant W(H). The asymptotic expression (87) includes this exponential 
growth, and also a (constant) pre-exponential factor. The pre-exponential factor, like W(H), is 
nonlocal: it is not determined by the behavior of the stochastic model in the vicinity of S and H. 
This is because the nominal 'frequency factor' K(H) can be computed only by integrating the 
system of ordinary differential equations (31), (32), (35), (36), (38) from S to H along the MPEP. 
The only exception to this is when 6' = —D^Wj^ for some differentiable function $, i.e., when 
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b is globally as well as locally gradient. In this case it is easily checked that W = 2$ and K = 1, 
so the pre-exponential factor in (87) is locally determined. 

The e ^ asymptotics which we have just derived for the exit location density p^ and the 
MFPT Er^ are quite familiar from the literature [45, 65, 69, 74]. In fact in the context of chemical 
physics the MFPT formula can be traced back to Eyring [29]. But our derivation of these formulae 
makes it clear that they are valid only when the classically forbidden wedge is absent. In general 
this will occur only in a single nongeneric case: when the drift field near H satisfies the local 
gradient condition, and moreover the eigenvalue ratio // > 1 (i.e., dih\H) < 0). In the generic 
case a wedge will be present, and the Hessian matrix Z{x) will not have a unique limit as x ^ H. 
As a consequence the inner approximation to the principal eigenfunction near H will have different 
Gaussian falloff rates as (x — H)/e^^^ ^ oo to either side of H, the asymptotic exit location density 
on do. will not be Gaussian, and equation (87) must at the very least be reinterpreted. 

Most previous work has concentrated on globally gradient models, for which b' = —D'Wj^. 
Such models are extremely nongeneric in that irrespective of //, they have no forbidden wedge: 
W = 2$ is always smooth at H, even when // > 1. It is perhaps this that has given rise to a 
general impression that the Hessian matrix Z(x) always has a well-defined limit as x ^ H, and 
that asymptotic exit location location distributions are always Gaussian. Generic models satisfy 
however neither b' = —D'Wj^ nor the local gradient condition c = 0, and we now see that their 
behavior is altogether different. In the next two sections we analyse the behavior of generic models 
with /i < 1 and // > 1 . 

7 Skewing and MFPT Asymptotics When d^b\H) > 

We now consider models with b\i(H) > 0, i.e., models in which the eigenvalue ratio // = 
\Xs(H)\/Xu(H) < 1. We shall show that generically, the asymptotic exit location distribution 
near H isa WeibuU distribution on the (9(e^/^) lengthscale, as mentioned in the Introduction. It is 
non-Gaussian and asymmetric. 

In this section we shall for simplicity take the linearized drift B'j(H) = b\j(H) to be 

/ Xu(H) \ / 1 \ 
B(H) = "' ' = I (88) 

V c X,(H) I \c -ijj 

as in Section 6, and take D(H) = I. The lower triangular form for B(H) can be arranged by an 
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appropriate linear change of coordinates near H, as can D(H) = I. Xu(H) = 1 can be arranged by 
a rescaling of time and noise strength. With this choice of B(H), the boundary dQ. will (near H) 
be parallel to the x^-axis. The most probable exit path (MPEP) from S to H will generically 
be tangent to dQ., and without loss of generality we assume that it approaches H from the first 
quadrant, as in Figure 2(a). 

The basic properties of // < 1 models were worked out in Sections 5 and 6. Generically 
there is a classically forbidden wedge emanating from H, on the boundary of which the Hessian 
matrix W^ij is discontinuous. If c < ('Subcase A') then within the wedge W^ij(x) has rank-1 
limit W^ij(H) as X ^ H, though if c > ('Subcase B') W(x) in the wedge is not expected 
to be a quadratic function of x near H. If x ^ H from outside the wedge, e.g., along the 
MPEP, then W^ij(x) -^ 0. Outside the wedge the action W is 'flat' near H. Since the outer 
approximation to the quasistationary density v^(x) contains an exponential factor exp {—W(x)/e), 
this is a sign that on the classically allowed side of H, the quasistationary density falls off only 
slowly (more slowly than quadratically). As a consequence the appropriate lengthscale for an inner 
approximation near H should be larger than (9(e^/^). On the (9(e^/^) lengthscale the asymptotics 
of u" and the exit location density p^ are those of (73)-(77), but an approximation on that lengthscale 
is not particularly useful. 

It was shown in Section 5 that when ji < 1, the tangent space T^h,o)M'^sq^ to the mani- 
fold M.'^sfi) C n X M^ at (i7, 0) generically equals A/'(+,+), the linear span of two eigenvectors of 
the linearized Hamiltonian flow at (i7, 0): the stable eigenvector (6^,0) and the unstable eigenvec- 
tor (e„,0). Since the MPEP q*, regarded as a trajectory in the phase space Cl x E^, terminates 
at (-ff, 0), it must approach (H^O) along the tangent vector (6^,0); since e^ lies in dQ., that is 
why it is tangent to dQ.. This tangency condition is of course an asymptotic statement, valid 
only as t ^ oo, i.e., as (H^O) is approached. In fact the 4-by-4 matrix T(H) representing the 
linearization of the Hamiltonian flow has both (6^,0) and (6^,^^) as stable eigenvectors. They 
have eigenvalues Xs(H) = —ji and —Xu(H) = — 1 respectively, so a more precise description of 
the t ^ oo asymptotics of the MPEP would be 

g*(t) — i7 ~ Cs exp(— /it) e^ -I- (7^ exp(— t) e^, t ^ oo. (89) 

Since // < 1, the first term is dominant as t -^ oo, and gives rise to the approach along dQ. 
However the coefficient C^, like Cg, is generically nonzero. The two coefficients can only be 
found by computing the MPEP explicitly: by integrating Hamilton's equations (31), (32) from 
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(«,/>) = (S^ 0) at t = — oo to (i7, 0) at t = oo. The fact that Cg is generically nonzero was taken 
into account when plotting Figures 2(a), 3(a), and 3(b); in those three figures, the slight deviation 
of the approaching MPEP from dO. is due to the e^ term. 

Explicit formulae for the eigenvectors e^ and e^ appear in (66)-(69); we may take e^ = (0, 1) 
and Bs = (fi — I, c). So (89) may be rewritten, if Aa; signifies x — H, as 

(Ax\ Ax^) ~ (Csip - l)e-\ CsC-^' + Csce-') . (90) 



We must have C's < and C's > 0, since it is our convention that the MPEP approaches H from 
the first quadrant. As t ^ oo and H is approached, the MPEP will generically be asymptotic to 
the curve Ax^ = A(Ax^Y, where 

A = Cs ( - — i 1 . (91) 

This asymptotic behavior occurs irrespective of the value of c; in this regard Subcases A and B are 
the same. 

Since the boundary layer near the characteristic boundary d^ will have thickness 0(e^/^), the 
inner approximation near H should be valid when Ax^ = 0(e^^^). But if the inner approximation 
is to be valid on a region containing a nontrivial (as e ^ 0) portion of the MPEP, it should also be 
valid when Ax^ = (9(e^/^). We see that x^ — H^ and x^ — H^ should be treated asymmetrically: 
when fj< I, the appropriate boundary region for the inner approximation is a strip near dQ. 
within which Ax^ ^ Ax^. This thin strip should extend from H along dO. in the direction of the 
approaching MPEP; equivalently, it should not extend in the direction of the forbidden wedge. 
As previously discussed, the appropriate lengthscale for an inner approximation on the wedge side 

Let us write (x, z) for (Ax\ (Ax^Y^'^] ; the change to noncovariant notation will emphasize the 
asymmetry. In the (x, 2;)-plane the boundary region will be the region where x, z = 0(e^^^). The 
MPEP will approach H = (0,0) as t ^ oo along an asymptotically linear trajectory z ~ A^^'^x. 
It is an easy exercise to show, using the matrix Riccati equation (38) in the linear approximation 
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near H, that irrespective of the asymptotic slope A^^^ the matrix of second derivatives 

Z{x, z) = 



dx^ dxdz 
V dzdx dz^ I 



(92) 



will have a finite (nonzero) limit as t ^ oo, z.e., as the MPEP approaches E. We accordingly 
expect that in the complement of the forbidden wedge, the action W near R behaves quadratically 
in X and z. 

This quadratic behavior will mandate Gaussian far-field asymptotics (as (x, z)/ e^^'^ -^ oo) for 
the inner approximation to u^, much as in the last section. The inner approximation is best written 
in terms of the stretched variables X = xj e^l^ and Z = zje^l^. It may be found by solving a 
linearized version of the approximate forward Kolmogorov equation £*/? = 0, as follows. In the 
linear approximation we take y (a;) sa B^ j(H)(x^ — H^) and D^^ (x) sa D'-' (if), and the Kolmogorov 
equation reduces to 

(e/2)D'\H)d,d,p - ddB',(H)(x^ - W)p] = 0. (93) 

Substituting both (88) and D(H) = I, and changing variables from (x\ x^) to X = (x^ — H^)/e^^^ 
and Y = (x^ - H^)/e^/^ yields 

In the e ^ limit this becomes 
and substituting Y = Z^ yields 

This partial differential equation for p = p(X^ Z) must be solved in the boundary region. 

In terms of the stretched variables (X, Z) we may view the boundary region as the right-half 
plane X > 0, and the boundary d^ as the Z-axis. The location of the classically forbidden wedge. 
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in terms of X and Z, is easily determined. As noted in Section 5 the wedge is bounded by the 
ray consisting of all multiples of e„; since e„ = (// + l,c) by (67), its boundary has equation 
Ax^ = [c/(p + l)]Ax\or Z = [c/(p + l)y/^e^^-^^/^^X^/^. In the e ^ limit the boundary of the 
wedge therefore becomes the X-axis: of the first and fourth quadrants in the (X, Z)-plane one is 
classically allowed, and the other is forbidden. We are taking the MPEP to approach (X, Z) = (0, 0) 
from the first quadrant, so it is the fourth quadrant that is forbidden. Equation (96) should therefore 
be solved only in the first quadrant. The extreme suppression of the quasistationary density in 
the forbidden wedge (on the (9(e^/^) lengthscale) allows us to set p = when Z < 0, for both 
Subcase A and Subcase B. 

A family of solutions of equation (96), each with far-field asymptotics that are Gaussian in X 
and Z, can be found by inspection. Each solution in the family is of the form 



p(X, Z) oc < 



Z^-i' exp(2BXZ - B^Z^), Z >0 

(97) 
0, Z <0 



for some constant B. Actually we shall use antisymmetrized (odd) versions of these solutions, since 
p must satisfy the Dirichlet boundary condition /j(0, •) = on account of absorption of probability 
on do.. Antisymmetrizing under x ^ —x yields 



p(X, Z) oc < 
Rewriting (98) in terms of Ax' = x' — W gives 



Z^-^ smh(2BXZ) exp(-B^Z^), Z >0 
0, z <0 



(98) 



'c(Ax2)(i/^)-i sinh [iBiAx^XAx^Y^^/e] exp [-B\Ax^)^^^/€] , Ax^ > 
0, Ax^ < 



(99) 



as the desired inner approximation for the case p < I, with C a constant to be found by matching 
to the outer approximation. Recall that in the far field, the inner approximation must match to 
an outer approximation of the form A"(a;)exp {—W(x)/e), where W is generated by incoming 
classical trajectories (including the MPEP). From (99) we can read off the behavior of W and K 
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as X ^ H when Ax^ > (i.e., in the complement of the forbidden wedge). Necessarily 

W(x) ~ -2B(Ax^)(Ax^y^^ + B\Ax^f^^ (100) 

K(x) ~ C(Ax2)(^/^)-^ (101) 

It is a useful exercise to verify that the formulae (100) and (101), irrespective of the choice of B, are 
consistent with the system of ordinary differential equations (31), (32), (35), (36), (38) in the linear 
approximation near H. The existence of a limiting value for the matrix Z(x, z) of (92) as x ^ H 
from outside the wedge is apparent: as a; ^ i7, 1^ is asymptotically quadratic in x and z. 

Notice that K ^ as x ^ H from outside the wedge: in particular, along the MPEP. In other 
words, the nominal frequency factor K(H) is generically zero when ji < \. This is an entirely 
new discovery, and has been confirmed in several models by numerical integration of the transport 
equation (36) from («,/>) = {S, 0) to (i7, 0). The fact that K(H) = is yet another reason why 
the classical formula (87) for the MFPT asymptotics cannot be generically applicable. 

It turns out that the constant B in the inner approximation is uniquely determined by the 
approach path taken by the MPEP. Since p, = d^W , differentiation of (100) yields 

p{H + Ax) ~ (-2B{Ax^yl^, -2Bijr\Ax^){Ax^)^^l^^-^ + lB^ijr\Axy^l^^-^) (102) 

near H, on the Ax^ = O ((Ax^Yj lengthscale. But in the linear approximation near H 

b(H + Ax)^(Ax\-jiAx^ + cAx^) (103) 

by (88), and substitution into Hamilton's equation (31) for x yields 

—(Ax\Ax^) ~ (-25(Ax2)i/^ + Axi,-25/i-\Axi)(Ax2)(i/^)-i + 25V"kAx2)(2/A')-i _^Ax2 + cAxi 

(104) 

as the equation of motion which must be followed by the classical trajectories near H giving rise 
to W. This equation is less complicated than it looks. When Ax^ = O ((Ax^Y), the dominant 
term in the second component on the right-hand side of (104) is the — // Ax^ term. The other terms 
may simply be dropped, and the equation simplifies to 

—(Ax\Ax^) ~ (-25(Ax2)^/^ + Ax\ -fiAx^) . (105) 

47 



It is trivial to verify that this asymptotic equation of motion near His compatible with the asymptotic 
approach path Ax^ = A(Aa;0^if andonly if i? = A~^^^. Since Ais the limitofthe ratio Aa;^/(Aa;^)^ 
along the MPEP as it approaches H, the constant B in the inner approximation may be computed 
numerically. 

Now that we have constructed an inner approximation to the quasistationary density v^ which 
is asymptotically accurate as e ^ 0, we may compute the asymptotic exit location distribution and 
MFPT asymptotics by the techniques sketched in Section 3. Since D(H) = I, equation (18) says 
that asymptotically, the density of the exit location measure on dO. is proportional to the normal 
derivative of u*'. This is simply the rate at which probability is absorbed on dQ., as a function of 
position. By differentiating the expression (99) with respect to Ax^ and setting Ax^ to zero we 
obtain (since B = A~^/^) 



i(- 



Pe(s) ~ < 



3(2/a*)-i exp[-(s/A)2/^/e], s > 0; 



0, s < 0. 



Here we have written s for Ax^ (the distance along dD., measured from H) for consistency with 
the Introduction and Section 6. The overall normalization of (106) is fixed by the condition that p^ 
have total unit mass. 

The asymptotic exit location density (106), which is localized on the s = (9(e^/^) lengthscale, 
is of an unusual form. It is the density of a WeibuU distribution [3], with shape parameter 2/// and 
scale parameter 1/Ae^/^. WeibuU-distributed random variables are simply powers of exponential 
random variables, and p^ may be viewed as the density of an 'offset' random variable 6^ equal 
to Af ^/^, where M is an exponential variate of mean A^/^e. The WeibuU distribution is decidedly 
'skewed'; in fact it is supported entirely on the s > side of the saddle point. That s > 0, rather 
than s < 0, appears here is solely a matter of convention. By convention the s > side of H is the 
side from which the MPEP approaches, as in Figure 2(a). For later use we note that 

E6, = / sp,(s)Js ~ Ar(l + /i/2)e^/^ e^O (107) 

Jo 

is the expected offset from H along d^, when // < 1, at the time of exit. 

The qualitative features of the asymptotic density (106) can be explained by reference to 
Figure 2(a). The quantity s^/^ in the exponent is roughly proportional to the square of the distance 
between the point (0, s) on dO. and the closest point on the approaching MPEP. This is consistent 
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with a picture developed elsewhere [56], according to which the MPEP is surrounded by a 'tube' 
of probability current, the tube having a Gaussian transverse profile. We have already discussed 
why the limiting pAs), on the (9(e^/^) lengthscale, is zero for s < 0. The points on dO. with s < 
are classically forbidden, and p^s) in the forbidden region falls to zero on the s = 0(e^^^) 
lengthscale, as summarized in (76), (77). On the (9(e^/^) lengthscale Subcase A and Subcase B 
differ significantly, but on the (9(e^/^) lengthscale they become identical. 

It is clear how the e ^ MFPT asymptotics may be computed from the inner approxima- 
tion (99). Asymptotically Xf\ i.e., (EtJ~\ is simply the flux of probability into dO.. As such it is 
proportional to the constant pref actor C in (99). If the inner approximation is to match up with the 
outer approximation, by (101) we must have 

C = ivexp(-W/(i7)/e), (108) 

where L is the limit of the quantity A"(a;)/(Aa;^)*^^/^^~^ as x ^ H along the MPEP. This formula 
for C has novel consequences. The factor exp {—W(x)/€) gives rise to the familiar Wentzell- 
Freidlin growth factor in Er^. But the pre-exponential factor in the asymptotics of (Er^)"^ is 
not proportional to the (generically zero) nominal frequency factor K{H), as in the classical 
formula (87), but rather to L, the rate at which K{x) approaches zero as the MPEP approaches H. 
In fact substituting the inner approximation into (22), and using B = A~^/^, yields 

(Er,)-^ ~ Af ~ ^^^i--^^detZ(5')exp i-W(H)/e) , e ^ 0, (109) 



which is the generic replacement for the classical formula (87) when // < 1. The A/det Z(S) 
factor arises from the denominator of (22), as in (87). We remind the reader that we are assuming 
Xu(H) = 1 and D(H) = I here. 

The generic applicability (when ji < 1) of this formula for the MFPT asymptotics, and the 
generic inapplicability of the classical formula (87), have not previously been recognized. It is 
remarkable that despite the nominal frequency factor K(H) equalling zero, the pre-exponential 
factor in (109) fails to be e-dependent. Naively one would have expected it to contain a positive 
power of e. A positive power of e is known to occur in the weak- noise reciprocal MFPT asymptotics 
of stochastic models where the nominal frequency factor equals zero on account of the exit location 
on do. converging to an unstable fixed point [54] . 
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8 Skewing When d^b\H) < 

We now consider models with b\i(H) < 0, i.e., models in which the eigenvalue ratio // = 
\Xs(H)\/Xu(H) > 1. The asymptotic exit location distribution near H is localized on the (9(e^/^) 
lengthscale, and as shown in Section 6 is generically non-Gaussian: it has different Gaussian falloff 
rates to either side of H. Although we shall not compute an explicit expression for its density, 
we shall work out a scheme for computing its moments of any desired order. Our treatment will be 
based on a probabilistic analysis, rather than on the construction of an inner approximation to the 
quasistationary density. 

In this section, as in Sections 6 and 7, we shall without loss of generality take the linearized 
drift5'j(i7) = b\j(H) to he 

and take D(H) = I. With this choice of B(H) the boundary dQ. near H will be parallel to 
the x^-axis, as in Figures 2 and 3. We also translate coordinates so that H = (0,0). With these 
normalizations the stochastic differential equation (1) becomes, in the linear approximation near H, 

dx\it) = x\it)dt + e^l^dwi{t) (111) 

dxlit) = -jix] dt + cxl(t) dt + e^/^dw2(t). (1 12) 

If the 'stretched' process (X(t), Y(t)) in the (X, y)-plane is defined to equal (^xl(t), xl(t)j /e^/^ 
we have 

dX(t) = Xdt + dwi(t) (113) 

dY(t) = -jiY dt + cXdt + dw2it). (114) 

We see that t i-^ X{t) is an inverted (repelling) Ornstein-Uhlenbeck process. On the (9(e^/^) 
lengthscale near H the region Cl becomes the right-half plane X > 0, and the fact that equation (113) 
does not involve Y indicates that with these normalizations the exit problem becomes essentially 
one-dimensional. 

Our interest is in the final approach to the boundary, which as e ^ will take place along 
the MPEP (most probable exit path) determined in Sections 5 and 6. Generically the MPEP, as 
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shown in Figure 2(b), is tangent to the stable ray e^ emanating from H. As computed in (68), 
Eg equals (// — 1, c). So as e ^ the final approach to dQ., in the linear approximation near H, 
will be increasingly concentrated on the line x'^ jx^ = c/(fj. — 1). On the (9(e^/^) lengthscale the 
approach path will not have a deterministic limit as e ^ 0. However the straight- line deterministic 
asymptotics should appear in the far field of the (9(e^/^) lengthscale; going backward in time from 
the boundary hitting time, the approach path should be asymptotic to the line Y/X = cjip — 1). 

The random variable 6^ = ^^^^Y{t,}, the displacement from H along dQ. at the boundary 
hitting time r^, is the quantity whose distribution we wish to compute. We define the time- 
reversed process (X(m), Y{u)\ m > 0, to equal iX{T^ — u), Yir^ — u))\ t^ is of course a random 
variable, which depends on the sample path. With this convention, X(0) = and YiO) = Q./e^l'^. 
A straightforward integration of (1 14) from u = Otou = T,i.e., from t = t, tot = t, — T, yields 
that for any T > 0, 

t-T 

Y(0) = Y(T)e-''^+ e-''''[cX(u)du-dw2(u)]. (115) 

JO 

As will shortly be seen, as e ^ the expected transit time of the final approach path tends to 
infinity. This justifies the taking of the T ^ oo limit when computing e ^ asymptotics. Taking 
the r ^ oo limit yields 

y(0)~ / e-^^[cX(u)du-dw2(u)], (116) 

Jo 

which is to be interpreted as a statement that the left and right-hand sides are distributed identically 
in the e ^ limit. But /q°° e~^^dw2(u)du is a Gaussian random variable of mean zero and 
variance 1/2//. It follows that 



6,/e^/2~ca(/i)-h3/V2/«, (iiv) 

where 3 is standard normal and the integral 



3{a)= / e-^"X(u)Ju (118) 

JO 

is a weighted area under the graph of the time-reversed process X{u), u > 0. The two terms 
in (1 17) are independent. 

The asymptotic exit location density p^(s) equals (d/ds) P {S^ < s}, the density of S^. In Sec- 
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tion 6 we deduced that when // > 1, p^is) has different Gaussian decay rates as s/e^^^ -^ ±00; 
from (78), 



Pe(s) ~ < 



exp 



/"(/" - 1) ,^2 



C2 + (/i - 1)2 

exp [-//(sVe)] , 



(sVe) 



s/e 
s/e 



1/2 



1/2 



+00, 



-00. 



(119) 



In deriving (119) the convention was adopted that the MPEP should approach H from the s > side; 
this amounts to assuming that c > 0. This being the case, it is interesting to compare (119) with 
the asymptotic equality in distribution (117). They are perfectly consistent: the Gaussian falloff 
of the density of the '5/^/2jI term in (117) may be viewed as the cause of the comparatively rapid 
Gaussian decay of pe('S) as s/e^/^ -^ —00, since the random variable 2f(/i) is non-negative. Infactby 
independence, the density p^(-) will be the convolution of the densities of e^/^2f(/i) and e^l^'il ^/Tjl. 
Equivalently, the generic asymptotic exit location density on the 0{e^l^) lengthscale will be the 
convolution of the density of c2f(/i) with a Gaussian (the density of'^|^/2jl). A striking conclusion 
follows: The skewing of the exit location distribution when fj> I is attributable to the asymmetry 
of the distribution of c3(fj.). Only if c = is this asymmetry absent. When c = 0, as we have 
already seen in Section 6, there is no skewing: the exit location distribution is asymptotically 
Gaussian, with variance 1/(2//). 

To determine the distribution of 3(fi), or at least its moments, we need to study the one- 
dimensional process t \-^ xl(t), conditioned on its exit at time r^. By the 'final approach path' we 
shall mean that segment of the trajectory t \-^ x^(t) which leaves some specified neighborhood 
of S and terminates on dQ. at time r^. Let a > be specified, and suppose that along the final 
approach path, xl (t) first reaches the point x^ = ae^l^aX time t = t^ — u; equivalently, that X{t) first 
reaches the point X = a at time t = t, — u. The process Xit), r^ — u < t < r^, is an inverted 
Ornstein-Uhlenbeck process conditioned to satisfy Xit) > for all t ^ {t,— u, t^), and X{t,} = 0. 

It is not difficult to compute the asymptotics of the distribution of u, the 'additional time to 
absorption' variable, in the large-a limit. The transition density p(Xo,to;Xi,ti) of an inverted 
Ornstein-Uhlenbeck process is of the form 



p(Xo,to;Xi,ti) = [27ra2_to]-^/'exp [-(X, - e"-'oXof/2al_,^] 



(120) 



where al = (e^^ — l)/2 is the variance at elapsed time z. If the process X(t) is conditioned to 
begin at a > at some specified time to, the probability of its having reached X = by time to + u 
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will by the method of images equal [16] 

2(1-/ p(a,to;X,,to + u)dXi}. (121) 

I JXi>0 J 

This absorption probability equals P {u < u}, the probability that the additional time to absorption 
is no greater than u. Substitution of (120) into (121), and some elementary manipulations, yield 
that 

P {u < u} ~ exp (_e-2(«-'°aa)^ ^ a ^ 00. (122) 

It follows that one may write 

u~loga + (5, a ^ 00, (123) 

where (5 satisfies 

P{(5 <u} = exp (-e-2") . (124) 



Equation (123) is an asymptotic equality in distribution, and the distribution of (5 is a so-called 
Gumbel (or double exponential) distribution of the sort that arises in extreme value theory [67] . 

It is noteworthy that if a is taken to equal Ce~^/^ for some C > 0, so that u is the amount of 
time that elapses between the moment the final approach path reaches x^ = C and the moment of 
final absorption at x^ = 0, the formula (123) implies that to leading order as e ^ 0, 

Eu ~ log (Ce-^/2) ~ (1/2) log(e-^) (125) 



irrespective of C We have commented elsewhere on the implications of this logarithmic growth [54, 
56]; see also Ludwig [50]. Its interpretation is as follows: the time needed for the process to make 
its final approach to the characteristic boundary grows logarithmically in e. This logarithmic 
growth is to be contrasted with the exponential growth of the MFPT Er^ as e ^ 0. The exponential 
timescale is the timescale on which a successful exit is expected to occur; when it occurs, however, 
it takes place on a much shorter (logarithmic) timescale. (1/2) log(e~0 is best viewed as the time 
needed for the deterministic MPEP to approach within an 0{e^l^) distance of the characteristic 
boundary. On that lengthscale (the width of the boundary layer) the MPEP ceases to be well- 
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defined; equivalently, the limiting approach path process ceases to be deterministic. The Gumbel 
random variable (5 is the additional amount of time needed for the process to reach the boundary. 
Here we are interested primarily in the implications of the asymptotic representation (123) for 
the time-reversed process X(u), u > 0, and the weighted area 2f(/i). The time-reversed process is 
an Ornstein-Uhlenbeck process satisfying 

dX(u) = —X(u) du + dw{u), (126) 

conditioned to satisfy X{0) = and X(u) > for all u > 0. Moreover, X(u) = a by definition. 
In the large-a limit we may write this condition as X(log a -i- (5) = a, so conditioning on the event 
(5 = M simply imposes an additional condition 

X(\oga + u) = a (127) 

onu \-^ X(u). In the language of random processes, once a value u for the random variable (5 is 
specified, the process X(u), u > 0, becomes a conditioned Ornstein-Uhlenbeck meander process. 
'Meander' refers to the fact that X(u) > for all u > 0, i.e., the fact that return to zero (i.e., to dO) 
is not allowed [23]. 

We shall shortly see that imposing the additional condition (127), and taking the a ^ oo limit, 
yields a well-defined process which we may denote Xu(u), u > 0. This being the case, define 

/•oo 

Mf^)= / e-^"Xs(u)Ju (128) 

to be the conditioned version of 3(fi). By (117) the moments of Gje^^^, the normalized displace- 
ment along do. at the time of hitting, satisfy 

E(&Je'^r ~ X:(^)^'(2/^r^'"'^'Ea(/.)'E3'=-' 

Jnn^-c^-iy^ik-iy.i] Ea(/i)' 
= E ( J c\2^i)-^^-'^l\{k -/)!!] y_^^^ Ea,(/i)' J[exp(-e-2«)] (129) 
where the integral arises from removing the conditioning (5 = u. This formula expresses the 
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sra-'^^^'- 



moments of the asymptotic exit location distribution p^(s) ds, when // > 1, in terms of those of 
the random variables 3u(fi). And by (128) the moments of 3u(fi) can be computed, by repeated 
integration, from the correlation functions (finite-dimensional distributions) of the process Xu(u), 
u>0. 

For any u, Xu(u), u > 0, is a Markov process whose transition probabilities may be computed 
by taking the above a ^ oo limit. However, it turns out to be time-inhomogeneous. It is preferable 
to express 3u(fj.) in terms of a closely related time-homogeneous process, which is based on 
Brownian rather than Ornstein-Uhlenbeck motion. This process is introduced as follows. Recall 
that a standard Ornstein-Uhlenbeck process o(u), u > satisfies 



o(u) = e-^w ((e^" - l)/2) (130) 



in the sense of equality in distribution; here w(t), t > 0, is a standard Wiener process. Formally, 
imposing the condition o(u) > for all u > is equivalent to imposing the condition w(t) > 
for all t > 0; in other words, the transformation (130) relates Ornstein-Uhlenbeck meander o"^ to 
Brownian meander w'^, just as it relates unconstrained Ornstein-Uhlenbeck motion o to Brownian 
motion w. Similarly, imposing the condition o'^(u) = a on Ornstein-Uhlenbeck meander o"^ 
at M = log a -I- M is by (130) equivalent to imposing the condition 



g-«.,.+ 



'w^ ((e^" - l)/2) = a (131) 

atu = log a -I- M on the corresponding Brownian meander w'^, i.e., imposing the condition 

w^ ([a^e^" - l]/2) = a^e". (132) 



By defining T = a^e^"/2 we see that for any fixed u, in the large- a limit this condition is to leading 
order a requirement that 

w^(T) = 2e-"r. (133) 

Informally, as a ^ oo and T ^ oo the conditioned Brownian meander corresponding to the 
original conditioned Ornstein-Uhlenbeck meander is forced by the condition (133) to drift upward 
at a mean speed 2e~". 

We define wUt), t > 0, to be the weak limit as T ^ oo of the standard Brownian meander 
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process w'^(t), t > 0, when constrained by the condition (133). We necessarily have, as an equality 
in distribution, 

l,(u) = e-"< ((e2« _ i)/2) , (134) 

so that 

/•oo 

Mp) = / e-^"Xs(u)Ju 
Jo 

= / e-(^+^^"u;t ((e2« - l)/2) Ju 

(l + 2t)-(^+^>/2<(t)Jt. (135) 



The last equality follows by a change of variables t = (e^" — l)/2. The formula (135) permits 
the computation of the moments of 3u(fi), as required by (129), from the correlation functions 
(i.e., finite-dimensional distributions) of the process w^. We shall shortly see that w^ is time- 
homogeneous, making this representation particularly useful. 

The n-point correlation functions of w^ may be computed by taking the T ^ oo limit of the 
n-point correlation functions of Brownian meander w'^, conditioned by (133). The evaluation of 
this limit is facilitated by the following fact. Recall that a three-dimensional Bessel process B(t), 
t >0, conventionally taken to satisfy B(0) = 0, is the radial coordinate in E^ of a diffusing particle. 
That is, B(t) equals [wl(t) + wl(t) + wl(t)y^^, where the Wi(t) are independent Wiener processes. 
It is a standard result [41] that the Bessel process B, when conditioned to satisfy B(t') = x' for any 
specified t' > and x' > 0, and Brownian meander w'^, when conditioned to satisfy w'^(t') = x', 
become identical in distribution on the time interval < t < t'. This allows us to substitute the 
Bessel process for Brownian meander, and to compute instead the T ^ oo limit of the n-point 
correlation functions of B, conditioned on B(T) = 2e~"r. The substitution of i? for w'^ simplifies 
the computation, for the Bessel process is (unlike Brownian meander) time-homogeneous. 

Denote by q(wi, ti;w2, ^2) the transition density of the standard Wiener process, i.e., 

q(wi,tuW2,t2) = [2Tr(t2-ti)r'^^ew[-i^2-w0V2(t2-ti)] . (136) 

If instants < ti < • • • < t„ are specified, the n-point correlation function p'^"^(-, ti; ... ; •, t„) for 
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the Bessel process B, which is defined by 



satisfies 



where 



p^'^\wi,ti; ... ■,Wn,tn) = P {B(t,) G w, + dw,, I < I <n} / Y[ dw. 



i=i 



n-l 



p^'''\wi,tu ... ; t«„, tn) = p^'^\wi,ti) Yl t(w„ tf, W,+l,t,+l) 



i=\ 



(137) 



(138) 



p^^\wi,ti) = J—^w^exp [-w\l2ti 



(139) 



and 



tiwi,ti;w2, t2) = (w2/wi)[q(wi, t^, W2, ti) - q(-wi,ti; W2, ti)] 



(140) 



are the probability density and transition density for the Bessel process. Conditioning on the event 
B(T) = 2e~"r, where T > tn, yields a process with n-point correlation function 



J^) 



p\'T{w]_,tu. . . \Wn,tn) 



/"^^Vi,ti;... ;u;^,t^;2e-"r,r) 
^(i)(2e-«r,r) 



(141) 



,(1) 



In particular, the conditioned density p)^ j{wi^ti) satisfies 



,(1) 



pt'iwut,) = p''Kw,,h) 



t(wi,ti;2e-^T,T) 
p(i)(2e-«r,r) 



(142) 



and the conditioned transition density tu^xiwi^tu W2, ^2) satisfies 



,(2) 



tu,T(Wi,tuW2,t2) 



PlTiwi,tuW2,t2) 



„(1) 



pVAw\,ti) 



t{Wi,ti\W2,t2) 



t(w2,t2;2e-^T,T) 
t(u;i,ti;2e-«r,r) 



(143) 



p;(l) 



,(1) 



het pr^'(wi^ti) and tu{wi^tuW2.,t2) be the T ^ 00 limits of p\j{wi^ti) and tu,Tiwi^ti\W2.,t2) 
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respectively. It follows by taking the T ^ oo limit of the two factors in brackets that 

P<."(,.„W = p">(„„.f,)e-«-/^(5i2!l=) 

V vwi J 



and 

iu(wuti;w2,t2) = e'" (*2-*'V2 / __^ \ [q(wi,tuW2,t2) - q(-wi,ti;w2,t2)]. (145) 

\smhvwiJ 

Here v = 2e~". p''l\w[,ti) and iu(wi,t[;w2,t2) are the density and transition density of the 
limiting process wl(t), t > 0. 

The transition density iu{wi^ti;w2.,t2) is invariant under time translation, so the limiting pro- 
cess is a time-homogeneous Markov process as promised. It is well known [62], and follows 
by examination of the transition density (140), that the three-dimensional Bessel process has gen- 
erator —(1/2)J^/Jt«^ — (l/w)d/dw. By examination of the formula (145) for t{i(t«i,ti; 1^2,^2)5 the 
process w^it), t > 0, has generator —(l/2)(f/dw^ — (v coth vw)d/dw. The function ucothutw is 
asymptotic tol/w asw -^ 0, and to u as t« -^ 00. This confirms that v = 2e~" has an interpretation 
as the strength of a superimposed upward drift. 

Now that the probability density and transition density of the process w^ are known, its n-point 
correlation functions p^^\wi,ti; ... ■,Wn,tn) follow from 

n-l 

p^^\wi,ti; ... ; t«„, tn) = p^i\wi,ti) Y[ h{w,, t,; w,+i,t,+i). (146) 

8=1 

Since 3u(fj.) is by (135) a weighted area under the process w^, its moments can be expressed 
in terms of these n-point correlation functions by repeated integration. It is not clear whether 
a closed-form expression for the distribution of 3u(p) exists. If it does it is likely to be quite 
intricate, as is suggested by the results of other authors. The problem of computing the distribution 
of the (unweighted) area under a Brownian bridge (i.e., a pinned Wiener process) was solved 
by Cifarelli [13] and Shepp [72], and the corresponding problem for a Brownian excursion by 
Louchard [47, 48]. More recently, Takacs [73] has computed the distribution of the integral of the 
absolute value of a Wiener process. All these distributions have closed-form expressions that are 
surprisingly complicated; they involve, for example, double Laplace transforms of the logarithmic 
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derivative of an Airy function. 

It is unfortunate that one cannot go directly from ( 144) and ( 145), or from the explicit expression 
for the generator, to a closed-form expression for the distribution of 2fs(/i). If such an expression 
were known, it would be possible to remove (by integration) the conditioning on (5 = u, and 
obtain a closed- form expression for the distribution of 3(fi). This in turn would yield a closed- form 
expression for the limiting exit location distribution on the (9(e^/^) lengthscale. But in the absence 
of such an expression one can at least compute the moments of the limiting distribution to any 
desired order, by using (129), (135), (144), and (145). 

The case of the first moment (the expected offset from H along dO., at the time dO. is reached) 
is particularly straightforward. By (129), 

E (6 Je^/2) ~ c E a(/i) (147) 

where 

E3(fi)= Eas(/i)J[exp(-e-2")]. (148) 

Ju=0 



Moreover, by (135) 

E 3dfi) = / (1 + 2t)-(^^3)/2 E wl(t) dt, (149) 

Jo 

in which 

Ewl(t)= / wp^l\w,t)dw, (150) 

Jo 

with the density p''i\w, t) given by (144). Evaluating the integral (150) yields 

E wld) = {v-' + vt)Qxf UvHJl] + J^e-^''^\ (151) 

which is a result of independent interest; this quantity is the expected distance from dQ. (on the 
(9(e^/^) lengthscale) at t time units before exit, if one conditions on the Gumbel random variable (5 
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equalling ii. Substitution of (151) into (149) and (148) yields, after various manipulations, 

E J(/i) = \ + -—=- (152) 

TT fi^ — 1 4y'Tr{fi — 1) 



where i?(-, •) is the Euler beta function. We conclude, by (147), that the expected offset from H 
at the time of exit, i.e., ES^ = f^ sp^(s)ds, is when jj, > I asymptotically equal to ce^/^ times this 
function of the eigenvalue ratio //. 

This result on the first moment provides additional information on the degree of skewing present 
in generic // > 1 models, over and above the differing sje^l'^ -^ ±oo asymptotics of ^^(s) given 
in (119). It is in fact possible to speculate, on the basis of the Gaussian falloff rates and the first 
moment, on the functional form of the generic asymptotic exit location distribution when // > 1. 
But we shall resist the temptation. 

We note briefly, in conclusion, that the probabilistic analysis of this section may be extended 
to models with // < 1 as well. If // < 1, it was shown in Section 5 that generically the MPEP is 
tangent to dO.. If B(H) is normalized as in (88), and D(H) = I, the MPEP approaches H = (0, 0) 
along a curve x^ ~ A(x^Y, where A can only be computed by integrating Hamilton's equations 
from S to H. In this case the conditioning on x^ = ae^/^ at time t = r^ — u, as e ^ 0, must 
be supplemented by a conditioning on x^ = A(ae^^^Y. It follows from the stochastic differential 
equation (112) that asymptotically, xI(t^) ~ xI(t^ — u)exp(— //u), irrespective of the value taken 
by the coefficient c. Since u ~ log a + (3, this implies that when // < 1, 

6, ~ A(ae^^yexp[-fi(\oga + (3)], (153) 

i.e., 

6, ~ Ae^/2exp(-/i0). (154) 

By examination, this offset random variable 6^ has the WeibuU distribution previously computed 
in (106) by the method of matched asymptotic expansions. For models with // < 1 , the probabilistic 
analysis and the method of asymptotic expansions are in agreement. 

Even though a probabilistic analysis can be performed when // < 1, the method of matched 
asymptotic expansions is superior in that it yields an asymptotic approximation to the quasi- 
stationary density in the boundary layer. Although our probabilistic results on the case // > 1 are 
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fairly strong, when // > 1 we have as yet no analogous boundary layer approximation. 

9 Two-Dimensional Ackerberg-O'Malley Resonance 

Now that we have to a large extent determined the generic asymptotic exit location distributions, 
we can examine the implications of our results for the singularly perturbed boundary value problem 
introduced in Section 2. Recall that if 

C, = -(e/2)D''d,dj - Vd, (155) 

is the generator of the process x^{t), t > 0, then the solution u^ of the boundary value problem 

C^u^ = in Q., u^ = f on dO. (156) 

will satisfy 

u,(y) = Eyf(x,(T,)) (157) 

where the subscript on E signifies the starting point ^^(O) for the process. As e ^ the function u^ 
is expected to 'level,' or tend to a constant, exponentially rapidly; this has been verified rigorously 
for the case of a non-characteristic boundary by Eizenberg [25]. The levelling is in agreement with 
the probabilistic picture that for any j/ G H, as e ^ it becomes overwhelmingly (exponentially) 
likely that a sample path for the process x^(t) will first flow toward S, and approach S to within 
an (9(e^/^) distance, before experiencing further fluctuations. As a consequence the expectation Ej, 
in (157) may up to exponentially small errors be replaced by E5, and u^(y) may be approximated 
by a J/ -independent constant. Since (if ^^(O) = S) the exit location on dQ. converges in probability 
to H, if / is continuous at H then u^(y) -^ f{H) for all y ^ Q.. 

Even though u^ will level exponentially rapidly and may be approximated by a constant, the 
constant itself will be e-dependent. It is possible to apply our results to determine its asymptotics, 
and the speed of its convergence to f(H), as e ^ 0. The key results are (107) and (152), in which 
we determined the expected offset from H along dO., at the time the process x^(t) exits Q.. In the 
standardization of the last two sections (the unstable eigenvalue Xu(H) taken to equal unity, and 
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D(H) taken to equal /), we found that the expected offset E&^ = J^ s p^i-s) ds has asymptotics 



E6, ~ [Ar(l + /i/2)]e 



A*/2 







(158) 



if /i < I, i.e., d.UiH) > 0, and 



E6e ~ c 



1 ^B(l/2,fi/2) 



TT fj/ 



1 40F(/i - 1) 



A/2 



(159) 



if /i > 1, z>., dib\H) < 0. Here // = \Xs(H)\/Xu(H) as always, A is a quantity that may be 
computed from the way in which the generically unique Wentzell-Freidlin trajectory from S to H 
(the MPEP) approaches H, and c is the off-diagonal element of the linearization ofbatH (see (88) 
and (1 10)). B(-, •) is the Euler beta function. 

It follows from (157), (158), and (159) that if / on dO. is continuously differentiable at H, then 
for all y ^ Q., u^(y) has leading e ^ asymptotics 



Ue(y) ~ < 



f{H)+[AY{l + p^ll)\nH)e^l\ 



f(H) + c 



TT fi^ 



i?(l/2,/i/2) 
1 ■ 4^(/i - 1) 



+ 



if /i < 1; 
f'(H)e'/\ if/i>l. 



(160) 



which are independent of y. If ji > 1 and c = (a local gradient condition at H, as discussed 
in Section 6) then the (9(e^/^) correction to f(H) will have zero coefficient. The same will occur, 
irrespective of //, if 6' = — i)'^$ ,,, i.e., if the drift field b is globally gradient. In these nongeneric 
cases the asymptotic exit location distribution will be a Gaussian centered on H, on the (9(e^/^) 
lengthscale, and the leading correction to f(H) in u^ will necessarily be o(e^/^). The same is 
true, incidentally, when D. has non-characteristic boundary. When Q. is attracted to S but dQ. is 
non-characteristic, it can be shown by the method of matched asymptotic expansions that the exit 
location distribution will generically be a Gaussian on the (9(e^/^) lengthscale, centered on the point 
on do. (generically unique) at which W attains its minimum. 

Theasymptoticsof(160)arestriking, especially in the case /i < 1. Since // need not be rational, 
the presence of a leading correction term proportional to e^/^ implies that u^ on Q. cannot in general 
be expanded in an asymptotic series in integral powers of e, or even in fractional powers. This 
has not previously been realized. We expect that a similar phenomenon will occur if the boundary 
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value problem (156) is generalized to read 

jC,u, = \'^QHn a, u, = f on da, (161) 

where Aq"^ (as defined in Section 2) is the limit as e ^ of the eigenvalue A^"^ of the n'th 
eigenfunction u" of C^. In general we should have u^ ~ C„(e)MQ, e ^ 0, where Mq is the limit 
of u" as e ^ 0, i.e., a solution of the degenerate first-order problem CqUq = Aq^^Mq. Since Aq ^ = 
and Uq is constant, this is a generalization from n = to arbitary n of the behavior we have just 
found. The small-e behavior of C„(e) for general n is likely to be at least as unusual as that of Co(e). 
To place these results in context, we remind the reader that the analogue for one-dimensional 
problems of the phenomenon we are investigating (the existence of a nontrivial e ^ limit on Ci. 
for the solution u^ of the singularly perturbed boundary value problem (156), or its generaliza- 
tion (161)), is known as Ackerberg-O'Malley resonance [1, 66]. In the one-dimensional case the 
partial differential equation C^u^ = Xu^ reduces to an ordinary differential equation, namely 

- (e/2)D(x)u'^ - b(x)u[ = Am,. (162) 

If an interval Q. = (xq, xi) is to have 'characteristic boundary' dO. = {xq, xi} in the sense of this 
paper, both xq and xi must be linearly unstable turning points; they must be zeroes of b(-) with 
b'(xo) and b'(xi) both positive. Moreover S G (xq, xi) must be a linearly stable turning point: the 
only zero of b(-) in Q., with b'(S) < 0. If A is arbitrary, irrespective of the choice of boundary values 
u^(x()) and u^(xi) the function u^ will normally converge to zero exponentially as e ^ 0, uniformly 
on any compact subset of Q.. The only exception to this occurs when A = Aq"^ for some n, i.e., 
when the boundary value problem is (asymptotically) 'in resonance' with the n'th eigenmode of the 
operator £,, equipped with Dirichlet boundary conditions. In this case u, converges to a nontrivial 
solution of the degenerate problem Cqu = Aq"^m. Since this convergence is uniform on compact 
subsets of (xo, xi), there are normally boundary layers near the endpoints xq and xi, of (9(e^/^) 
width as e ^ 0. 

The explicit solution of the n = one-dimensional problem is instructive. In one dimension 
the drift field b(-) is necessarily a gradient; b(x) = —D(x)dx^(x), where 

^(x)= r[-b(y)/D(y)]dy, (163) 

Js 

and the action W equals 24> as usual. When n = we set A = Aq*^ = 0, and the differential 
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equation (162) simplifies. The closed-form solution of (162) is (cf. O'Malley [66]) 



u^(x) 



/:■ exp {W(y)/e} dy 
III exp {W(y)/e} dy 



Ue(xQ) + 



Jxo 



exp {W(y)/e) dy 



rxi 
I xo 



exp {W(y)/€) dy 



u,(x,). (164) 



If W(xo) < W(xi) (resp. W(xi) < W(xo)), then by examination u^(x) -^ u^(xo) (resp. u^(x) -^ 
u^(xi)), uniformly on compact subsets of Q as e ^ 0. Moreover, u^ both levels and converges 
exponentially rapidly. The exponential levelling, and the constant limit, have the same interpretation 
in terms of the stochastic exit problem as they do in two-dimensional models. 

The one-dimensional boundary value problem cannot usually be solved in closed form when 
n > 1, but a method of matched asymptotic expansions [66], involving the construction of both 
inner and outer approximations to u^, may be employed to show that for any x G H, u^(x) has 
an asymptotic expansion in fractional powers of e. Although this is more complicated than the 
asymptotic behavior for the one-dimensional n = problem, it is still simpler than the behavior 
we have discovered in its two-dimensional analogue. In two-dimensional Ackerberg- O'Malley 
resonance (at least, when Q. has characteristic boundary) the outer expansion of u^ in the body of ^ 
may involve non-fractional powers. 

The presence of irrational powers of e in the outer expansion of u^, as in (160), suggests that 
they may also be present in the outer expansions of the eigenfunctions u" and u" of C^ and C*. 
For this reason, unlike many authors we have refrained from approximating the quasistationary 
density u*' in the body of Q by a formal asymptotic series, since it is unclear what powers of e 
should be present. Instead, we have worked only to leading order. As we noted at the beginning 
of Section 4, for an outer expansion to be useful it must match to an inner expansion. And the 
expansion beyond leading order of the quasistationary density v^ in the boundary layer, if dO. is 
characteristic, remains an unsolved problem. 

10 Conclusions 

The generic features of the two-dimensional stochastic exit problem, when exit from the region Q. 
occurs near a saddle point H, are now clear. As the noise strength e ^ 0, the exit location near H 
will be concentrated on the (9(e^/^) lengthscale near H (if the eigenvalue ratio // < 1) or the 0(e^/^) 
lengthscale (if // > 1). In the // < 1 case the exit location distribution is asymptotic to the WeibuU 
distribution (106), which includes a scale factor that can only be computed from the approach path 
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taken by the MPEP (the optimal, or most probable trajectory) from S to H. In the // > 1 case the 
limiting exit location distribution, whose moments are computable (see, e.g., (152)), contains no 
free parameters: it is determined by the behavior of the stochastic model in the vicinity of H. 

In both cases the limiting distribution will be 'skewed': non-Gaussian and asymmetric. Nor- 
mally, it is Gaussian only when the deterministic drift b satisfies 6' = —D^Wj^ for some potential 
function $, or when // > 1 and a local version of this equality (the c = condition of Section 6) 
holds near H. These cases, which are characterized by the absence of a classically forbidden 
'wedge' emanating from H, are nongeneric. Although our two-dimensional stochastic model 
differs from the barrier crossing models employed in chemical physics, we believe that the generic 
skewing phenomenon is related to the phenomenon of 'saddle point avoidance' [2, 4, 46] . A number 
of authors have in fact already noted the presence (in particular models) of a classically forbidden 
region. In the literature the boundary of the forbidden region is sometimes called the 'stochastic 
separatrix,' or a 'switching line' [2, 4, 24, 40, 65]. 

It is clear from our treatment that the generic features of models with // < 1 are particularly 
interesting. In such models the nominal frequency factor K(H) (the value of the WKB prefactor 
atx = H, which would normally be interpreted as a factor by which the frequency of excursions 
to the vicinity of H is multiplied) equals zero. This feature, like the anomalously large lengthscale 
over which the exit location distribution is spread [(9(e^/^) rather than (9(e^/^)] can be traced to the 
unusual approach path taken by the MPEP. When // < 1 the MPEP is generically tangent to the 
separatrix dQ. at H, as in Figure 2(a). This grazing behavior causes the exit location distribution to 
be anomalously wide. It also causes the WKB prefactor K(x) to tend to zero as x ^ H along the 
MPEP, as can be shown by integrating the system of ordinary differential equations (31), (32), (35), 
(36), (38) from S to H. Another unusual feature of generic models with // < 1 is that the classical 
formula (87) for the mean exit time asymptotics must be replaced by (109). The formula (109) 
is unaffected by K(H) equalling zero, and by the fact that generically, the Hessian matrix of the 
Wentzell-Freidlin action W does not exist at H. 

The system (31), (32), (35), (36), (38), which generates the WKB (outer) approximation to 
the quasistationary density, is particularly useful for numerical work. We commented briefly on 
the interpretation of the matrix Riccati equation (38) for didjW in terms of symplectic geometry. 
It turns out, as we may explain elsewhere, that the transport equation (36) for K(-) also has a 
geometric interpretation. The natural setting for the outer approximation K(x) exp {—W(x)/€) is 
the theory of semiclassical expansions for partial differential equations [22, 59], which has deep 
geometric underpinnings. 
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It should be possible to extend the approach of this paper in several different directions. It has 
recently become clear that caustics (folds in the unstable manifold A^"^ q^ in phase space comprising 
the most probable fluctuational trajectories, as in Section 4) occur very frequently [24, 55]. Their 
effect on first passage and exit time phenomena is now under investigation. Another extension, 
of particular value in applications, would be to the case of singular diffusion. The results of 
this paper apply to what is known in physics as the overdamped limit of barrier crossing models; 
the analysis of exit location distributions in models with underdamped dynamics will require an 
extension to the case when the diffusion tensor is allowed to become singular [65]. An extension 
to higher dimensionality would also prove useful, and is now under study. 

As regards the connections with the phenomenon of Ackerberg-O'Malley resonance discussed 
in Section 9, we cannot resist quoting O'Malley [66]: 

We note that resonance is sometimes related to certain exit-time problems for 
stochastic equations, but regret to report that the mathematical phenomenon under 
discussion (despite much attention in the literature) has not yet substantially helped us 
understand much new physics. 

It now appears that the desired connections to physical models can indeed be found, by going from 
one-dimensional resonance (for ordinary differential equations) to multidimensional resonance 
(for partial differential equations). Skewed exit location distributions and saddle point avoidance 
have implications for the asymptotic expansions used in analysing multidimensional resonance, 
and vice versa. 
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